{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 21\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### With air resistance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next we'll add air resistance using the [drag equation](https://en.wikipedia.org/wiki/Drag_equation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I'll start by getting the units we'll need from Pint."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "kilogram"
      ],
      "text/latex": [
       "$\\mathrm{kilogram}$"
      ],
      "text/plain": [
       "<Unit('kilogram')>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I'll create a `Params` object to contain the quantities we need.  Using a Params object is convenient for grouping the system parameters in a way that's easy to read (and double-check)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>height</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.0025 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.019 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>18.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "height                      381 meter\n",
       "v_init             0.0 meter / second\n",
       "g             9.8 meter / second ** 2\n",
       "mass                  0.0025 kilogram\n",
       "diameter                  0.019 meter\n",
       "rho         1.2 kilogram / meter ** 3\n",
       "v_term            18.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "params = Params(height = 381 * m,\n",
    "                v_init = 0 * m / s,\n",
    "                g = 9.8 * m/s**2,\n",
    "                mass = 2.5e-3 * kg,\n",
    "                diameter = 19e-3 * m,\n",
    "                rho = 1.2 * kg/m**3,\n",
    "                v_term = 18 * m / s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can pass the `Params` object `make_system` which computes some additional parameters and defines `init`.\n",
    "\n",
    "`make_system` uses the given radius to compute `area` and the given `v_term` to compute the drag coefficient `C_d`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Makes a System object for the given conditions.\n",
    "    \n",
    "    params: Params object\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    diameter, mass = params.diameter, params.mass\n",
    "    g, rho = params.g, params.rho, \n",
    "    v_init, v_term = params.v_init, params.v_term\n",
    "    height = params.height\n",
    "    \n",
    "    area = np.pi * (diameter/2)**2\n",
    "    C_d = 2 * mass * g / (rho * area * v_term**2)\n",
    "    init = State(y=height, v=v_init)\n",
    "    t_end = 30 * s\n",
    "    dt = t_end / 100\n",
    "    \n",
    "    return System(params, area=area, C_d=C_d, \n",
    "                  init=init, t_end=t_end, dt=dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's make a `System`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>height</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.0025 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.019 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>18.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>0.0002835287369864788 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.4445009981135434 dimensionless</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>y             381 meter\n",
       "v    0.0 meter / secon...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>30 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>0.3 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "height                                              381 meter\n",
       "v_init                                     0.0 meter / second\n",
       "g                                     9.8 meter / second ** 2\n",
       "mass                                          0.0025 kilogram\n",
       "diameter                                          0.019 meter\n",
       "rho                                 1.2 kilogram / meter ** 3\n",
       "v_term                                    18.0 meter / second\n",
       "area                         0.0002835287369864788 meter ** 2\n",
       "C_d                          0.4445009981135434 dimensionless\n",
       "init        y             381 meter\n",
       "v    0.0 meter / secon...\n",
       "t_end                                               30 second\n",
       "dt                                                 0.3 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the slope function, including acceleration due to gravity and drag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Compute derivatives of the state.\n",
    "    \n",
    "    state: position, velocity\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: derivatives of y and v\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    rho, C_d, g = system.rho, system.C_d, system.g\n",
    "    area, mass = system.area, system.mass\n",
    "    \n",
    "    f_drag = rho * v**2 * C_d * area / 2\n",
    "    a_drag = f_drag / mass\n",
    "    \n",
    "    dydt = v\n",
    "    dvdt = -g + a_drag\n",
    "    \n",
    "    return dydt, dvdt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As always, let's test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0 <Unit('meter / second')>, -9.8 <Unit('meter / second ** 2')>)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the same event function as in the previous chapter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Return the height of the penny above the sidewalk.\n",
    "    \"\"\"\n",
    "    y, v = state\n",
    "    return y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then run the simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here are the results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0.0</th>\n",
       "      <td>381 meter</td>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.3</th>\n",
       "      <td>380.559 meter</td>\n",
       "      <td>-2.913855777777778 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.6</th>\n",
       "      <td>379.2553998560887 meter</td>\n",
       "      <td>-5.676321799192868 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>0.9</th>\n",
       "      <td>377.15535917269835 meter</td>\n",
       "      <td>-8.166374221525693 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1.2</th>\n",
       "      <td>374.35521895425103 meter</td>\n",
       "      <td>-10.311720070623485 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                            y                                   v\n",
       "0.0                 381 meter                  0.0 meter / second\n",
       "0.3             380.559 meter   -2.913855777777778 meter / second\n",
       "0.6   379.2553998560887 meter   -5.676321799192868 meter / second\n",
       "0.9  377.15535917269835 meter   -8.166374221525693 meter / second\n",
       "1.2  374.35521895425103 meter  -10.311720070623485 meter / second"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>y</th>\n",
       "      <th>v</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>21.300000</th>\n",
       "      <td>20.486949941310563 meter</td>\n",
       "      <td>-17.99999999508715 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21.600000</th>\n",
       "      <td>15.086949942543686 meter</td>\n",
       "      <td>-17.99999999642989 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21.900000</th>\n",
       "      <td>9.686949943439785 meter</td>\n",
       "      <td>-17.99999999740564 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22.200000</th>\n",
       "      <td>4.286949944090969 meter</td>\n",
       "      <td>-17.999999998114706 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22.438164</th>\n",
       "      <td>0.0 meter</td>\n",
       "      <td>-17.99999999852377 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                  y                                   v\n",
       "21.300000  20.486949941310563 meter   -17.99999999508715 meter / second\n",
       "21.600000  15.086949942543686 meter   -17.99999999642989 meter / second\n",
       "21.900000   9.686949943439785 meter   -17.99999999740564 meter / second\n",
       "22.200000   4.286949944090969 meter  -17.999999998114706 meter / second\n",
       "22.438164                 0.0 meter   -17.99999999852377 meter / second"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "results.tail()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final height is close to 0, as expected.\n",
    "\n",
    "Interestingly, the final velocity is not exactly terminal velocity, which suggests that there are some numerical errors.\n",
    "\n",
    "We can get the flight time from `results`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "22.438163885803732 second"
      ],
      "text/latex": [
       "$22.438163885803732\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "22.438163885803732 <Unit('second')>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "t_sidewalk = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the plot of position as a function of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saving figure to file figs/chap21-fig01.pdf\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3hUVf7H8fekN3oLHaQceq8iEIKuddVFwLUhiqtCKCLlh2V1WQsuHQQElbWuZe29EnqRGlrg0EsCofcQSsjvjzvokIUYMMmdJJ/X8+TJzL13ki9xzCfn3nPP15ORkYGIiIi/CXC7ABERkYtRQImIiF9SQImIiF9SQImIiF9SQImIiF9SQImIiF8KyutvaIwpDqwCnrHWvmmMCQEmAV2BdGCstXaEz/HdgReB8sBsoKe1dm9e1y0iInnLjRHUVKCiz/PhgAFqAC2B+40xPQCMMfWA6UBPoBSwEfggL4sVERF35GlAGWPuB4oCq3023w+8YK09ZK3dBowGHvHuuxf4ylo7z1qbBjwBtDPG1MrDskVExAV5dorPGFMdeBa4Gvjeu604zqm7RJ9D1wMNvY/rAUvP77DWphpjdnr3b8zm9w3FGZntxjmFKCIi/iMQJweWWGtP+e7Ik4AyxgQC7wKDrbUpxpjzu6K8n1N9Dk8FInz2++7LvD87WgJzL6tgERHJa+2Beb4b8moE9XfAWms/zbT9hPdzuM+2COC4z/5wLuS7Pzt2A/znP/8hOjr6Ml4mIiK5LSUlhXvuuQe8v6t95VVA/RWoYIzp4n1eBJgCtAJScCZJJHv31eG3U36J3n0AGGMigCpceErw96QDREdHU6lSpSutX0REctf/XILJk4Cy1tbxfW6MSQDGe6eZHweeNcaswjmlNxiY4D30PWCeMSYGWAiMAFZYazfkRd0iIuIef7hR9xlgDbAWWAJ8gjMVHWvtauBB7/P9QH2gmztliohIXsrzG3UBrLVNfB6nAXHej4sd+wlOaImISCHiDyMoERGR/6GAEhERv6SAyoaMjAy3SxARKXRcuQaVn7z73To+mrGBkOBAIsKCCA8NIiIsmOhSkVQqG0XlskWoVC6KSmWjCA4KdLtcEZECQwH1O4KDnUFm2ul00k6nA85KHBt3Hr7guJCgAOpUK0nDmqVpWKM0tauUIDhIA1QRkSulgPodd15r6BZbm7TTZzl56iypaWc5lnqa3ftPkLT3ODv3HCNp7zGS951g1ab9rNq0H4CwkEBa1YvmmiYVaF6nHCHBGl2JiFwOBVQ2BAR4iAgLJiIsmFLFnG31qpe64Jgjx0+xZssBVntDaueeY8xJSGZOQjLhoUG0bhBNbPPKNK5VhoAAjwv/ChGR/EUBlUOKRYXSrlEF2jWqAMCeg6nMS0hm7spkNicdYdayJGYtS6J86UhubFuNzi2rUDQyxOWqRUT8lwIql5QrGcEdsbW4I7YWu/YdZ/aKZH78ZTu795/g31+t5Z3v1tG+SUW6xNSkavmibpcrIuJ3FFB5oEKZKO76k6F751osXbeHbxduY4XdS/zSncQv3UmretF0ja1F3eol3S5VRMRvKKDyUGBgAK0blKd1g/Ls3n+CL+Zs5qdftrM4MYXFiSnUq16Sv15naFK7DB6PrlOJSOGmedAuKV86kke7NGL603+i+7W1iQwPJnHrQZ55dSFPTJnPms373S5RRMRVCiiXFS8Syn031uXfT19Hj5vqEhUezNotB3hiynz+Pm0BG3YccrtEERFX6BSfn4gIC6Zb59rcdHV1vpizmc9nbyZhwz4SNuyjfZOK9LipLtGlIt0uU0Qkz2gE5Wciw4O5+/o6vP7UddzRqSYhQQHMTUim979m8NoXqzl64rTbJYqI5AkFlJ8qGhlCz1vq88qwzsS2qEz6uQy+nLOFh1/8ic9nb+bM2XNulygikqsUUH6ubIkIBt7VjPEDY2hSuwwn0s4y/cs19Bs9k6Xr9rhdnohIrlFA5RNXVSzGPx9uyzO9WlOhdCTJ+44z/PVF/OO1hSTtPeZ2eSIiOU6TJPIRj8dDy3rRNKldlm/mb+H9Hy3L1u9l5caZ3NahBndeZwgP1X9SESkYNILKh4KDAri9Y02mDbuW61pV4Wx6Bp/M3ETvf81g7opkNVgUkQJBAZWPFS8SSv87mzJmQAdqVi7OgSNpjHx3KU9PXcDOPTrtJyL5mwKqAKhdpQSj+3egb7fGFIkIYdWm/fQfM5O3v00k7fRZt8sTEbkiCqgCIjDAw/VtqjF1WGf+1LoqZ9Mz+GjGRuJGzWRxYorb5YmIXDYFVAFTNDKEft2bMLJve6qVL8reg6k8N/0XXnxzMQeOnHS7PBGRbFNAFVB1q5dk/MCO9Lq1AeGhgSxcvZve/4rny7mbST+nSRQi4v8UUAVYYGAAt3esweQhnWnTIJqTp87y2udrGDxxDpuSDrtdnohIlhRQhUCZEuE89UBrnnqgFaWLh7Np52EGjZ/N9C/XcPKUJlGIiH9SQBUibRqUZ8rQWG7rUAOAz2dvpu+oeC2ZJCJ+SQFVyISHBvHQbQ0YM6AjV1Usxt5DJxn++iL+9fYSDh1Nc7s8EZFfKaAKqZqVizN2QAd63Vqf0JBA5q3cRe+R8fz0y3atRCEifkEBVYg5kyhqMnlILM3rlOXEyTNM/G8CT09dwK59x90uT0QKOQWUUK5kBM8+1IZB9zSnWJSzEkXf0TP5aMYGzqar75SIuEMBJYCzUnpMs0pMGeo0SDxz9hxvf7uOx8fPZuPOQ26XJyKFkAJKLlA0MoSBdzXjuUfaUq5kBFt3HWXwhDlM/3INaZqSLiJ5SAElF9WkdlkmDe7E7R19pqSPnsnKDftcrkxECgsFlFxSWGgQvW5twOgBHaheoSh7Dqby9LQFTPhgBcdTT7tdnogUcAoo+V21Kpdg7GMdue/GugQHBfDzkh30HhnP/JW73C5NRAowBZRkS1BgAN2vrc3EQTHUq16Sw8dO8dLbS3jxzcUc1A2+IpILFFByWSqVLcKIPtfwaJdGhIcGsXD1bvroBl8RyQUKKLlsAQEebm5XnclDYmlRt9yvN/g+M20hKQdOuF2eiBQQCii5YmVKhPNMr9YMuqc5RSJCSNi4j76jZ/LFHPWcEpE/TgElf8j5G3xf+b9YOjStyKnT6bz+xRr+7+W5bE856nZ5IpKPKaAkRxSLCmXIvS34+4OtKVUsDLvjEI+NncX7P1rOnNVySSJy+YLy8psZY24BXgSqA3uBkdbaacaYEGAS0BVIB8Zaa0f4vK6793XlgdlAT2vt3rysXbKnVf1o6l9Vije/SeT7hdt474f1LFi1i37dm1C7Sgm3yxORfCTPRlDGmPLAx8D/WWuLAN2A8caYZsBwwAA1gJbA/caYHt7X1QOmAz2BUsBG4IO8qlsuX2R4MHFdG/NC76spXyqSbbuPMmTiHN74ai1pp7VckohkT54FlLV2N1DGWvudMSYAJ2zOAseA+4EXrLWHrLXbgNHAI96X3gt8Za2dZ61NA54A2hljauVV7XJlGtUsw8TBMfwlpiYAn87aRP8xs1i9eb/LlYlIfpCn16CstceMMRHAKeBHYDKwD+fUXaLPoeuBht7H9Xz3WWtTgZ0++8WPhYUE8eCf6zOqfweqRhdh9/4TPDllPlM+Xklq2hm3yxMRP+bGJIk0IBLnVN6DwADv9lSfY1KBCO/jqEz7Mu+XfKB2lRKMGxjD3X8yBAV6+G7hNuJGzWTpuj1ulyYifirPA8pae85ae9pauxR4FWjh3RXuc1gEcL6l64lM+zLvl3wiOCiAu66vw7iBMdSqXJz9h08y/PVFjHlvGUdPaPFZEblQXk6S6GiMWZZpcyhwCEjBmSRxXh1+O62X6LvPe4qwCheeEpR8pFr5oozq154HbqlPSFAAs5Yl0WfkDOYmJGu5JBH5VV5OM08AKhpjHgcmAK2BXsBfcALqWWPMKpxTeoO9xwC8B8wzxsQAC4ERwApr7YY8rF1yWGBgAF061aRNw2he/m8CazYfYOQ7S5mzIpredzSmZNEwt0sUEZfl5Sy+I8BNQBfgIM7pvYestbOBZ4A1wFpgCfAJMNX7utU416qmAvuB+jhT1KUAqFA6ihcebUefro0JDw1i0ZoU+oyM5+fFWnxWpLDzFPRfAsaYasDWGTNmUKlSJbfLkSzsO3SSKZ+s/HXiRNPaZejbrQllS2o+jEhBlZSUROfOnQGqe28z+pWWOhK/cX7x2cfvbkaRiGBWbNhH3Kh4vpm3hXNafFak0FFAiV/xeDx0al6ZyUNjade4Ammn05n62WqemDKP5H2auClSmCigxC+VKBLGsB4tebJnS0oUCSVx60H6jZ7Jx/EbSU/X4rMihYECSvxa24YVmDI0ls4tK3Pm7Dne+iaRwRPnsHXXEbdLE5FcpoASvxcVEcJjf23G8L+1pUyJcDYlHWHguNm8+/06zpxNd7s8EcklCijJN5rVKcukwZ24uV110s9l8OFPG3hs3Gw27DjkdmkikgsUUJKvRIQF82iXRozo047ypSPZkXKMIRPnMP3LNWrlIVLAKKAkX2pQozQvD+5EF28rj89nb1YrD5ECRgEl+VZocCAPqJWHSIGlgJJ873wrj7vUykOkQFFASYEQHBTA3d5WHjV9WnmMe385x1LVykMkP1JASYFSrXxRRvdrzwO31CMkKID4pTvpMzKeBat2uV2aiFwmBZQUOE4rj1pMHNyJ+leV4vCxU4x4awkvvbWEQ8fS3C5PRLJJASUFVsUyUbzYux2PdmlEWEgg81ftIm5kPDOX7VQrD5F8QAElBVpAgIeb21Vn8pBYmtYuw7HUM4x9bzn/nP4L+w6ddLs8EclCtjvqGmPKAc2BskA6Thfc5dbaA7lUm0iOKVsyguEPt2XGkh28/uValq7bQ9yoeB74c32ub12VgACP2yWKSCZZBpQxJgi4G3gMaAycBg4BgUBJ7zG/AFOAD6y1WmZa/JbH4+HaVlVpVqccr3yykkVrUpjy8UrmrkimX/cmlC8d6XaJIuLjkqf4jDEdgVVAD2A6UBuIsNZWsNaWA0KApsB7QF9gvTEmJtcrFvmDShYN48merRh6XwuKRYWwevN++o6eyeezN5GuxogifiOrEdQg4E5r7eqL7bTWZgBrvB9TjDFNgX8Cs3K6SJGc5vF4aN+kIo1qlub1L9Ywa3kS079cy7yEXfS7swlVo4u6XaJIoecp6LOZjDHVgK0zZsygUqVKbpcjfmpxonO678CRNIICA/jrdbW5I7YWQYGaRySSm5KSkujcuTNAdWvtNt99lzNJIgKoDoRm3metXf4HaxRxVat60dQfUoo3vl7LD4u28+7361mwajf972xCjUrF3S5PpFDK1p+Hxph7gb0416SWZvpYkmvVieShyPBg+nZrwvOPXE25khFs2XWExyfM4e1vEzl9Ro0RRfJads9fjMCZKHEVUD7TR4XcKU3EHY1rl2HS4E78uf1VZGRk8NGMjQwYO4v12w66XZpIoZLdU3xFgUnW2u25WYyIvwgLDeLh2xtyTeMKTPwwgaS9xxk6aS5/vuYq7ruxLmGh2T47LiJXKLsjqHeAnrlYh4hfqle9FBMHxdCtcy08Hg9fzt1CvzEzWblxn9uliRR42f0zcBSw3BhzD7ANuOCGXGttbA7XJeI3QoID6XFTPa5uVIGJH65g666jPD11Ade3qcoDt9QnMjzY7RJFCqTsBtQ7wHHgGyA198oR8V81KxVn7GMd+WTmRj74cQM/LNrO0nV76NO1Ma3qRbtdnkiBk92Aagm0ttauys1iRPxdUGAAd15raNOgPBM/XMGGHYd5bvovxDSrxN9ub0jRyBC3SxQpMLJ7DcoCuhlExKtqdFFG9utAr1vrExIcyKzlScSNjGfeymS18hDJIdkdQY0A3jTGTAI2A2d8d1prv83pwkT8XWCAh9s71qRV/Whe/m8CazYf4F9vL6Vtw/I82qURJYuGuV2iSL6W3YB63/t59EX2ZeCsbi5SKFUoHcULj7bjh0XbeOPrRBau3s2qTft56NYGdG5ZGY9HrTxErkS2AspaqwXJRLIQEODhxqur06JuNJM/TmDZ+r1M+HAFcxOSievWmLIlItwuUSTf+b12G5fFGKPp5lKolSkRzrMPtWHgXc2ICg9mud1L31HxfDN/K+fUykPksmQ1ghpojBkGTAR+ttaeudhB3qaGt+D0hEoF4nO8SpF8xOPxENuiMk1rl2HqZ6tYsGo3Uz9dxdyEZPp3b0KFMlFulyiSL2TZbsMY8xfgH0BVnD5Pa4H9gAcog9Nlty2wA3jOWvtx7pZ7+dRuQ9w2f+Uupn66isPHTxESFMA9N9Tlto41CFSbeZErb7dhrf0M+MzbKfcmnDAqh7OSRAqwDBhhrZ2b82WLFAztGlegYc3SvPbFamYtS+KNr9cyf1Uy/bs3pWp5NUYUuZTsTpKYhTrlilyxopEhDLq7OR2bVmLyRwls2HGYx8bNovu1hq6xtQgO0jwkkcz0f4VIHmpRtxyTh8ZyQ9tqnE3P4L0f1vP4+Nls2nnY7dJE/I4CSiSPRYQFE9e1MS/0vproUhFs232UQRPn8ObXa9UYUcSHAkrEJY1qluHlQZ24vWMNMjIy+GTmJvqPmUXi1gNulybiFxRQIi4KCw2i160NGNmvPZXLRZG87zjDJs9j2merOHnqrNvlibgq221BjTFlgUZAMM40819pLT6RP6ZO1ZJMeDyGD37awMfxG/l63lYWJ+6hX7fGNKld1u3yRFyRrYAyxvQCpuCEU2Zai08kBwQHBXLfjXVp16gCEz5cwZbkI/x92kKua1WFB29tQJQaI0ohk90R1BDgNeAJa+2xK/1mxpjrgJeAWsBeYJS1dpoxJgSYBHQF0oGx1toRPq/rDrwIlAdmAz2ttXuvtA4Rf3ZVxWKMGdCBT2du4v0fLT8t3sGy9Xvpc0cjWjco73Z5Inkmu9egKgMT/mA4VQY+AZ7H6S11FzDCGHM9MBwwQA2c5oj3G2N6eF9XD5gO9ARKARuBD660DpH8ICgwgO7X1mbioBjqVC3BwaNpPP/GYka9u5Qjx0+5XZ5InshuQP0IdP6D36sa8J619jNr7Tlr7RKcm3/bAfcDL1hrD3mXuhgNPOJ93b3AV9baedbaNOAJoJ0xptYfrEfE71UuV4SX+rbnodsaEBIcyJwVyfQZGc/cFWqMKAVfdk/xrQTGGmNuBTYAp313WmuH/t4X8C6H9OuSSMaYkkB74B2cU3eJPoevBxp6H9cDlvp8nVRjzE7v/o3ZrF8k3woM8HBbhxq0qhfNpI8SWLVpPyPfXcrsFdH0vqMRpYqFu12iSK7IbkB1BH4BwnEWiPV12X/GGWOKAV96v+Yy7+ZUn0NSgfMNdKIy7cu8X6RQKF86kucfvZofFm3n31+t5Ze1KazZcoCHbq1P55ZV1BhRCpzsrsXXKae+oTGmNvAFzojpHpzQw+czOOFz3Pv4RKZ9mfeLFBoej4cb2lZzlkz6eCVL1+1hwocJzFmRTFy3JpQrqb/bpOC4nPugyuH0fKqPc+1qHfCatXbLZXyNDjjhNBV40lqbAaQZY1JwJkkkew+tw2+n/BK9+85/jQigCheeEhQpVEoXD+eZXq2ZvTyJVz9fzYoN++g7Kp77b67HTVdXJ0CtPKQAyNYkCWNMK5xrT3/B6Qe1D6dJ4SpjTItsfo0awNfAM9baJ7zhdN47wLPGmNLe/k2DvdsA3gNuM8bEGGNCgRHACmvthux8X5GCyuPxENO8MpOHxtKucQXSTqcz7bPVPDFlHsn7dIJB8r/sjqDGAO8DvX2DxRgzCRgFZOcUYBxQBGdq+Qif7ZOBZ7zfYy1OaL6KM8rCWrvaGPOg93lFnOtW3bJZt0iBV6JIGMN6tGTBql288ukqErcepP/omdx9fR1u71iDwECtaCb5U5Yddc8zxpwEmlhrbabtBlhmrfXbHtbqqCuFybHU07z+xRril+4EoGbl4gy4synV1BhR/FRWHXWz+6fVbpz7mDK7Crjim3dFJGcViQhh4F3NePahNpQuHs6mnYcZOG4W7/2wnjNnz7ldnshlye4pvneAV40xjwGLvNvaAuP47VqRiPiJFnXLMXlIJ978OpHvFm7j/R8tC1fvpv+dTahVuYTb5YlkS3ZHUC/grCbxXyAJZ7bd+8BHwFO5U5qI/BERYcH06dqYF3u3o3ypSLbtPsrgCXN446u1nFJjRMkHshVQ1trT1tq/AaVxRk6NgeLW2sHW2jO5WaCI/DENa5Zm4uAYbu9YA4BPZ22i/+iZrN2ixoji3y55is8YcxPwk7X2jPdxZpWdORLqByXi78JCnMaI7ZtUZMKHK9iRcoxhk+dxc7vq9LipLhFhauUh/iera1BfA9E4bTG+zuI49YMSySdqVynB+IEd+fDnDXw8YyPfzN/KksQU4ro1oZlRY0TxL5cMKGttwMUei0j+FhwUyL03/NYYcXPSEZ59dSHXtqxCr1vrExUR4naJIkD2V5KIN8YUv8j2MsaYZRd7jYj4t+oVijGmfwfuv7kewUEB/LxkB3Gj4lm0ZrfbpYkAWV+DisFpdQHOauaPGGMy3/NUF6fJoIjkQ4GBAXSNrUXr+tG8/N8E1m07yAtvLKZ9k4o88peGFIsKdbtEKcSyugZ1AGdNPI/3Iw6nHft5GTgrig/KtepEJE9ULleEEXHX8M38Lbz97TrmJiSzcuM+Hr69IR2aVlQrD3FFVtegVuOsFIExZibQxVp7KK8KE5G8FRjg4db2vzVGXLlxP6P/s8zp4ttVjREl72V1ii/CWnu+UeDN57dd7Fif40Qkn4suFclzj1zNj7/s4N9frWFxYgprR+7ngT834E+t1RhR8k5WkySOGfPrvNPjOGvuZf44v11EChCPx8P1baoyZWgsrepFcyLtLJM+SuCZaQtJOXDC7fKkkMjqGlQscND7OMc66opI/lGqWDhPP9iKOSuSmfbZahI27qPv6Jn0uKkuN7e7ikA1RpRclNU1qNkXewxgjAkBGgEbrLVHc688EXGbx+OhY7NKNK5Vhlc/X83chGRe+3wN8xJ20a97EyqXK+J2iVJAZfc+qJrGmNnGmDbe61CLvR/bjTFtcrVCEfELxYuEMvS+FjzZsxUlioSybttBBoydxUczNpCerlYekvOyu0LEyzjXmrYB9wGVAAO8AozNlcpExC+1bVieKUNj6dyyMmfOnuPtb9cxaOIctu464nZpUsBkN6DaAwOttSnA7cA31tqNwGtAk9wqTkT8U1RECI/9tRnD/9aWMiXC2Zx0hIHjZvPu9+s4c1atPCRnZDeg0oBgY0wkzqoS33m3RwP6s0mkkGpWpyyTBnfi5nbVST+XwYc/beCxcbPZsEO3TMofl92A+gFntPQJkAp8ZYzp7N32ZS7VJiL5QERYMI92acSIPu2oUDqSHSnHGDJxDtO/XEPa6bNulyf5WHYD6hFgKc5I6mZr7QmgJTALeCx3ShOR/KRBjdJMHNyJLjE1Afh89mb6j5nFms37Xa5M8itPRkbGZb3AGFMUCLDWHs6dknKWMaYasHXGjBlUqlTJ7XJECoUNOw4x8cMVbE9x7uO/8epq9Ly5nhojyv9ISkqic+fOANWttdt892W7z5MxprcxZidwCDhgjNltjBmWo5WKSIFQu0oJxg2M4e4/GQIDPHy3YBtxo2aybP0et0uTfCS790ENBl7CmW7eHugAjAOGGmMG5F55IpJfBQcFcNf1dRg3sCM1KxVj/+GT/OO1RYx7fznHUk+7XZ7kA1ktdeQrDnjUWvu+z7b5xpjtwPPAhByvTEQKhOoVijG6fwc+n72Z//ywnvilO1lu99Lnjka0bVjB7fLEj2X3FF8ZYMlFti/DuWlXROSSAgMDuCO2Fi8P7kS96iU5fOwUL765hJfeXsKhY2lulyd+KrsBtQbodpHtdwLrc64cESnIKpaJYkSfa3jkLw0JCwlk/spdxI2MZ9aynVzuhC0p+LJ7iu8Z4BtjTFtgoXdbW+AGoEtuFCYiBVNAgIdbrrmKlt7GiAkb9jHmveXMSUgmrmtjNUaUX2VrBGWt/RHoDJzCWYuvK3AUaGmt/Tr3yhORgqpcyQj++XBb+ndvQmRYEEsS99BnZDw/LNqm0ZQA2R9BYa2dA8zJxVpEpJDxeDxc17oqzeqU5ZVPVvHL2hQmfbSSOSuS6de9CdGlIt0uUVyUZct3YDzOaOkU8BkwTP2fRCSnlSoWzlMPtGJewi6mfraKVZv203f0TO67sS63XKPGiIVVVqf4hgN/BkbitNS4GWftPRGRHOfxeGjftCJThsbSoWlFTp1O5/Uv1jBs0lx27jnmdnnigqwCqitwt7X2JWvtKJxZfLcZY7RWiYjkmmJRoQy5twVPP9CKkkVDWb/9EP3HzOK/P2/grBojFipZBVQlLpxCvsR7fLlcrUhEBGjdoDyTh3bmulZVOJt+jne+W8egCXPYkqwOP4VFVgEVCPzaecxam4FzLSokt4sSEQGICg+m/51N+efDbSlbIpwtyUd4fPxs3vlOjRELg2wvFisi4pampiyThsRyi7cx4n9/3sCAsbOx2w+6XZrkot+bZt7TGHM80/H3GmMuaPBirZ2S45WJiPgIDw3ikS6NuKZJRV7+7wp27jnGkJfncmv7Gtx7Yx3CQrJ914zkE1n9F90B9M60LQV4INO2DEABJSJ5ov5VpZgwqBPv/7Cez2Zt4os5m1m8NoV+3ZvQsGZpt8uTHHTJgLLWVsvDOkREsi00OJCet9SnXeMKTPwwgW27j/LkK/O5oW01HrhFjRELCl2DEpF8q1blEox9rCN3X1+HoEAP3y/cRtzIeJauU2PEgkABJSL5WnBQAHf9yTB+YAy1Khdn/5E0hr++iLHvLePoCTVGzM8UUCJSIFQtX5RR/TvwwC31CQkKYOayJOJGxjN/1S63S5MrpIASkQIjMMBDl041eXlwJ+pfVYrDx0/x0ltLGPHWYjVGzIdcmZdpjGkFfG2tLet9HgJMwlleKR0Ya60d4XN8d+BFoDwwG+hprd2b54WLSL5QoUwUL/Zux3cLt/HWN2tZsGo3qzft56HbGtKpeSU8Hi0+mx/k6QjKGOMxxjwE/PLM/PAAABB9SURBVMiFK1IMBwxQA2gJ3G+M6eF9TT1gOtATKAVsBD7Iw7JFJB8KCPBwc7vqTBocS9PaZTiWeoZx7y9n+OuL2HfopNvlSTbk9Sm+4Tj3Vj2fafv9wAvW2kPW2m3AaOAR7757ga+stfOstWnAE0A7Y0ytPKpZRPKxsiUjGP5wWwbc2ZTI8GCWrd9L3Kh4vlu4jXPn1BjRn+V1QE211jYHlp7fYIwpjnPqLtHnuPVAQ+/jer77rLWpwE6f/SIiWfJ4PFzbqgpThsbStmF5Tp46y5SPV/L01AXs3n/C7fLkEvI0oKy1F5tOE+X9nOqzLRWI8NmfyoV894uIZEvJomE8cX9L/q9HC4pFhbB6s9MY8fPZm0nXaMrv+MMsvvN/voT7bIsAjvvsD+dCvvtFRLLN4/FwTeOKTB4SS0yzSpw+k870L9fwf5PmsiNFDcP9iesBZa09hLPGn/HZXIffTusl+u7ztqKvwoWnBEVELkuxqFAG3dOcv/dqTcmiYdjthxgwdjYf/mTVGNFPuB5QXu8AzxpjShtjqgGDvdsA3sPp5BtjjAkFRgArrLUb3ClVRAqSVvWimTI0luvbVOVs+jne/X49g8bPYXPSYbdLK/T8JaCeAdYAa3E6934CTAWw1q4GHvQ+3w/Ux2k/LyKSIyLDg+nbrQnPP3I15UpGsGXXER6fMIe3v03k9Bk1RnSLJyOjYF8Y9I7Its6YMYNKlSq5XY6I+Lm0U2d557t1fDVvCxkZUKlsFAPubEqdaiXdLq1ASkpKonPnzgDVvbcZ/cpfRlAiIn4hLDSIv93ekJfirqFimSiS9h5n6KS5vPbFatJOnXW7vEJFASUichH1qpdi4qAYunWuhcfj4cs5W+g7eiYrN+5zu7RCQwElInIJIcGB9LipHmMGdKB6haLsOZjK01MXMOmjBE6cPON2eQWeAkpE5HfUrFScsY915N4b6xAUGMAPi7YTNyqexYkpbpdWoCmgRESyISgwgDuvNYx/vCO1qxTnwJE0npv+C2P+s4wjx0+5XV6BpIASEbkMVaOLMrJfB3rdWp+Q4EBmLU8iblQ881YmU9BnRec1BZSIyGUKDPBwe8eavDw4hgY1SnHk+Gn+9fZSRry1hINH1RgxpyigRESuUIXSUbzwaDv63NGI8NAgFq7eTZ+R8cxYskOjqRyggBIR+QMCAjzceHV1Jg3pRLM6ZTlx8gzjP1jBP15fxN5DmRsxyOVQQImI5ICyJSL4x0NtGHhXU6LCg1m+fi99R8Xz7YKtaox4hRRQIiI5xOPxENvCaYx4daPynDyVziufrOKpqfPZtV8dgi6XAkpEJIeVKBrGE/e3YliPlhSPCmXN5gP0Gz2Lz2ZtUmPEy6CAEhHJJe0aV2Dy0FhimjuNEf/91VqGvjyH7WqMmC0KKBGRXFQ0MoRBdzfn2YfaULpYGBt2HOaxsbN4/0fLmbNqjJgVBZSISB5oUbcck4bEckPbapxNz+C9H9bz+PjZbNqpxoiXooASEckjkeHBxHVtzAu9rya6VATbdh9l0MQ5vPn1WjVGvAgFlIhIHmtUswwvD+rErR2uIiMjg09mbqL/mFkkbj3gdml+RQElIuKCsNAg/nZbQ0b2bU/lclEk7zvOsMnzePXz1ZxUY0RAASUi4qo61UoyfmAM3a+tjcfj4au5TmPEhA173S7NdQooERGXhQQHct+NdRk7oANXVSjG3oOp/H3aQiZ+uILjhbgxogJKRMRP1KhUnDGPdeC+G+sSFBjAT4t3EDcynsVrC2djRAWUiIgfCQoMoPu1tZk4KIY6VUtw8Ggaz/37F0a9u7TQNUZUQImI+KHK5YrwUt/2PHRbA0KCA5mzIpk+I+OZu6LwNEZUQImI+KnAAA+3dajBpMGdaFSzNEdPnGbku0t54Y3FhaIxogJKRMTPlS8dyfOPXk1c18aEhwbxy9oU+oyM5+fF2wv0aEoBJSKSD3g8Hm5oW43JQ2JpUbccJ06eYcKHCTz76kL2HiyYjREVUCIi+UiZEuE806s1j9/djCIRwazYsI+4UfF8M29LgWuMqIASEclnPB4PnZpXZvLQWNo1rkDa6XSmfraaJ1+ZT/K+gtMYUQElIpJPlSgSxrAeLXni/pYULxLK2i0H6D96Jp/O3Eh6ev5v5aGAEhHJ565uVIEpQ2OJbVGZ02fP8cbXiQx+eS7bdufvxogKKBGRAqBIRAgD72rmNEYsHs6mnYcZOG4W7/2wPt82RlRAiYgUIC3qlmPykE7c6G2M+P6PlsfHz2bjzkNul3bZFFAiIgVMRFgwfbo25sXe7ShfKpJtu48yeMIc3vhqLafyUWNEBZSISAHVsGZpJg6O4faONQD4dNYm+o+eydot+aMxogJKRKQACwsJotetDRjZrz2VyxVh1/4TDJs8j2mfrvL7xogKKBGRQsBULcmExzty57W1CQzw8PX8rfQdFc8K67+NERVQIiKFRHBQIPfeWJdxAztSo1Ix9h46yTOvLmTCBys4nnra7fL+hwJKRKSQqV6hGGP6d6DHTXUJDgrg5yU7iBsVz6I1u90u7QIKKBGRQigwMIBunWsz4fEY6lYrycGjp3jhjcWMfMd/GiMqoERECrHK5YowIu4a/nZ7A0JDApmbkEzvf8Uze3mS6608FFAiIoVcYICHW9s7jREb1yrNsdTTjP7PMp7/92IOHDnpWl0KKBERASC6VCTPPXI1fbs1ISIsiMWJKcSNjOeHRe40RlRAiYjIrzweD9e3qcrkIbG0rFeOE2lnmfRRAs9MW0jKgRN5Wku+CChjTGNjzEJjzAljzGpjTEu3axIRKchKFw/n7w+2ZtA9zSkSEULCxn30Gz2Tr+bmXWNEvw8oY0wI8AXwIVAceAH40RhT1NXCREQKOI/HQ0yzSkwZGkv7JhVJO53Oq5+vZtjkeSTtPZbr39/vAwqIAYKtteOttWestR8Aa4E73S1LRKRwKF4klKH3teDJnq0oUSSUddsO0n/MLD6asSFXGyPmh4CqB6zLtG090NCFWkRECq22DcszZWgsnVtW5szZc7z97TqGTZ7HmbO5s0J6fgioKCA107ZUIMKFWkRECrWoiBAe+2szhv+tLWVKhLNl11FS03Jn0dmgXPmqOesEEJ5pWwRw3IVaREQEaFanLNOGXcvJU2cpGhmSK98jP4ygEgGTaVsd73YREXFJcFBAroUT5I8R1EzAY4wZCEwC7gAaAZ+5WpWIiOQqvx9BWWtPAzfiBNNB4CngdmvtPlcLExGRXJUfRlBYa9cA17hdh4iI5B2/H0GJiEjhpIASERG/pIASERG/lC+uQf1BgQApKSlu1yEiIpn4/G4OzLyvMARUeYB77rnH7TpEROTSygObfTcUhoBaArQHdgO5s2CUiIhcqUCccFqSeYfH7Z7zIiIiF6NJEiIi4pcUUCIi4pcUUCIi4pcUUCIi4pcUUCIi4pcUUCIi4pcUUCIi4pcUUCIi4pcKw0oSV8wY0xiYitPBdwvwoLX2f+52LsyMMQ8C04BTPpvjrLVvuVSS3zDGtAK+ttaW9T4PwekK3RVnVZOx1toRLpboFy7ycwoFjgGnfQ5bYK39kxv1uckYcx3wElAL2AuMstZOKyzvJQXUJXjfAF8A44EOOB19fzTGVLXWHnW1OP/SDBhjrR3mdiH+whjjAXoBozPtGg4YoAZQDPjeGJNsrX07j0v0C1n8nBoCB6210Xlflf8wxlQGPgHux/ld1Bz4wRizDYihELyXdIrv0mKAYGvteGvtGWvtB8Ba4E53y/I7zYEEt4vwM8OB3sDzmbbfD7xgrT1krd2G84v5kTyuzZ9c6uek95SjGvCetfYza+0579mbWUA7Csl7SQF1afWAdZm2rcf5604AY0wgzunP+4wxu4wxm4wxw7x/GRdmU621zYGl5zcYY4rjLIiZ6HNcYX8//c/PyasZUNYYs8oYs8cY85ExpqIL9bnKWjvXWvvo+efGmJI4C1+voJC8lxRQlxYFpGbalgpEuFCLvyqD88vlLaA6zvnw3t6PQstau+sim6O8n33fU4X6/XSJnxPACWA+0BnnNNZJ4LO8qssfGWOKAV8CvwDLvJsL/HtJ16Au7QQQnmlbBHDchVr8krU2BejosynBGPMyzvW6Ke5U5bdOeD/7vqf0froIa+3jvs+NMY8D+4wxla21O10qyzXGmNo416ASgXv47T1U4N9LGkFdWiLOX2++6nDhsLpQM8bUN8YMz7Q5BEhzox5/Zq09BKRw4XtK76eLMMb80xhT12dTiPdzoXtfGWM64IyaPge6WmvTCtN7SSOoS5sJeIwxA3Gmc96Bc72lUJ9qyOQwMMgYkwRMB5oC/YG+rlblv94BnjXGrMI55TcYmOBuSX6pEdDCGHO39/kE4Btr7T4Xa8pzxpgawNfAU9balzPtLhTvJY2gLsFaexq4ESeYDgJPAbcXtv9JsmKtTQZuxZk9dBRnSuxz1tqPXS3Mfz0DrMGZDboE5+c11dWK/FMv4BCwCdiGcz/UfW4W5JI4oAgwwhhz3OfjXxSS95I66oqIiF/SCEpERPySAkpERPySAkpERPySAkpERPySAkpERPySAkpERPySbtQVySHGmDdxVpm+lOE4q1HPBIpYa/NkaRrvor7zgR7W2g1ZHBcALALus9bavKhNJCsaQYnknAE4q0yXx2nXAtDKZ9toYIH38YmLvD639AdWZhVOANbac8A/KYA3fEr+pBt1RXKBMaYBsBqo7u3X41YdYcAOINZauyabr9kM9LLWzsrN2kR+j07xieQhY0wMPqf4jDEZwF3AEziLfy4F7gWG4CzvcxR4wlr7jvf1RYAxOK1NMoB4YEAWrSv+Chz2DSdjzN+Bh3HapawDnrTWfufzms9wRoOzcuCfLHLFdIpPxH0vAY8BbYAqwHKcYGoJfApMM8ac7yf1Kk6QXY/T6iQDpw34pf7YvBn4/vwTY8xfvN/rXpwVsL8BPjLGFPV5zffAtVl8TZE8oYAScd9ka+1Ma20CzurVx3FGNRYYi9P3p7ox5iqcEdHd1tol3lHRfTitwW+4xNdugbOg6HnVgFPAdu+px38CXYAzPsck4qyQXSdH/nUiV0h/IYm4b5PP41Rgm7X2/MXh8z2QQoGq3sfWmAtalUXgjKq+vsjXLgfs93n+Ls5Mwy3GmGU4XVrfsNae9DnmgPdz2cv8d4jkKI2gRNx3JtPzc5c4Lsh7bFOgic9HbeCNS7zmHOA5/8TbLqY5zohrAdATWOWd1HHe+d8L6dn+F4jkAgWUSP6xDggGIq21m6y1m4DdwCickLqYFJzJEAAYY7oAj1hrf7TWDsAZeR0DbvJ5TRmf14q4Rqf4RPIJa601xnwJvG2MiQP2AS/gTK5Yf4mXLQMa+zwPBEYZY/bgzBhsA0R7H5/XmN8aBoq4RiMokfzlfpww+Rynk2ox4Dpr7eFLHP8Nzmw/AKy1HwHP4oy6NgDPA32ttfE+r+kAfG+t1Sk+cZVu1BUpwIwxETht02+w1i7PxvEBwHacmYJzc7k8kSxpBCVSgFlrU3FGS3HZfMltwBaFk/gDBZRIwTcOaGQyzU3PzDt6egp4NE+qEvkdOsUnIiJ+SSMoERHxSwooERHxSwooERHxSwooERHxSwooERHxS/8P1tFu/4aIzOUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_position(results):\n",
    "    plot(results.y)\n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Position (m)')\n",
    "    \n",
    "plot_position(results)\n",
    "savefig('figs/chap21-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And velocity as a function of time:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZxcVZn/8c+trl7T2UlIwpoE8gAKkbAIEVnFEZElgApGRgZcYHRExPEXxVEUoyjooOIMuAwiLlEUlN0AAuJKQBIIgUcSTCBk35fu9FJVvz9udVLpVHVXd1fVrer6vl+venXfU6duPZ1XpZ8+5z73nCCVSiEiIlJuYlEHICIiko0SlIiIlCUlKBERKUtKUCIiUpbiUQdQTsysHjgGWAkkIg5HRKQa1ADjgXnu3pb5hBLU7o4Bnow6CBGRKvRW4I+ZDUpQu1sJ8NOf/pRx48ZFHYuIyKC3atUqZs6cCenfv5mUoHaXABg3bhz77rtv1LGIiFSTPS6rqEhCRETKkhKUiIiUJSUoEREpS0pQIiJSliqmSMLMpgK3AEcArwCXuvu8LP32B34IHAesAf7D3R8oZawiIjJwFTGCMrM64LfAL4ARwGxgrpkNy9J9DvAcMBr4EDDHzCaVKlYRESmMikhQwMlArbvf5O4d7j4HeAF4b2YnM5sCHA183t3b3f33wD3AZaUIMtnWSvv610vxViIig16lJKjDgBe7tb0EHJ6l36vuvr2XfkWx5p5vsfzWT9CxaU0p3k5EZFCrlATVDLR0a2sBmvrZryhSiQSkkrSvXlqKtxMRGdQqJUFtBxq7tTUB2/rZryhqR4wFoHOzRlAiIgNVKQlqEWDd2g5Jt3fvt7+ZNfbSryji6QSlKT4RkYGrlDLzx4DAzK4CbgbOJyw3vzuzk7u7mS0AZpvZZ4DpwDnA8aUIMj5cIygRkUKpiBGUu7cDZxAmpg3ANcC57r7WzGaaWeYU3vnAoYT3QP0AuMzdF5Yizp1TfBpBicggd/XVVzN79uydx4lEgunTp/O3v/2tYO9RKSMo0knmhCztPwV+mnH8GmEyK7muEVTH5rWkUimCIIgiDBEZRFbOmU3rkr+X7P0aJ09j/IXX9NpvxowZzJo1i1mzZlFTU8Of/vQnGhoaOPbYYwsWS0WMoCpFrLGZoK6RVFsLyR0lqcsQEYnE9OnTicViO0dM9957L2eddVZB/zCvmBFUJQiCgNoRY2hf8yqdm9ZS0zg06pBEpMLlM5qJQiwW4+yzz+bee+/lyCOP5JFHHuGuu+4q7HsU9GxCfPjeAHRsXh1xJCIixTVjxgweeeQR5s6dy5QpU5g4cWJBz68EVWDxEWMAFUqIyOA3efJkDjjgAG666SbOOeecgp9fCarAakeEIyglKBGpBjNmzGDt2rW8853vLPi5laAKLD48HEHpZl0RqQYzZ85k4cKFjBgxouDnVoIqsHjXCEo364qIDIgSVIHVpkdQnel7oUREpH+UoAos1jCEWEMzqY42ki1bog5HRKRiKUEVwa5FY1VqLiLSX0pQRRDPmOYTEZH+UYIqgl2l5hpBiYj0lxJUEajUXERk4JSgiqBWpeYiIgOmBFUEWu5IRGTglKCKYPd9oZIRRyMiUpmUoIogVtdArGkYJDpJbN0YdTgiIhVJCapIdl2HUqm5iEh/KEEVya5KPpWai4j0R0XsqGtmVwJXAqMBB6529ydz9P0SMAtoz2h+l7s/Xuw4M3WtJqERlIhI/5R9gjKz84BPA6cDLwEfAO4zs4PcPdtv/2nAx939lhKGuYfadKGEbtYVEemfSpjiGw98xd0XuXvS3W8DEsDhOfofBcwvWXQ57FyPTyMoEZF+KYsRlJnVAaOyPJVy9+9263si0Ay8kOU844FxwCwzOx5YD9yQTmoltXOKTyMoEZF+KZcR1HRgZZbH65mdzOyNwC+Az7l7tt/8Y4EngJuBfYErgJvM7MzihZ7drgVj15FKJkr99iIiFa8sRlDpAoagpz5m9i7gDuB6d/96jvMsAE7OaHrCzO4AzgPuL0iweYrF66hpHkli20Y6t67feU1KRETyUy4jqB6lq/h+DnzQ3b/WQ78TzOwT3ZrrgB3FjC+XXdN8WvJIRKSvymIE1RMzew/wFeBUd/9bL91bgevN7GXgQeBU4H3pryVXO3wsbcs9TFAHRBGBiEjlKvsERXhPUz3wqJlltl/o7veZ2S3AAe5+hrs/Y2YXA18Hfgm8Clzi7k+VPGognl5NQjfrioj0XdknKHef1svzl3c7vhO4s6hB5al2ZNfGhZriExHpq4q4BlWpdo6gNmoEJSLSV0pQRVQ7chwAnZtWRRyJiEjlUYIqopqhI6EmTmL7ZpLtrVGHIyJSUZSgiigIYtSq1FxEpF+UoIosPiKc5uvYqGk+EZG+UIIqsq5KPpWai4j0jRJUkXVV8nWqkk9EpE+UoIqsq5JPpeYiIn2jBFVktV0jKE3xiYj0iRJUke3cuHDTGm27ISLSB0pQRRara6BmyAhIdpLYuiHqcEREKoYSVAnEVcknItJnSlAlULtzTT7dCyUiki8lqBKId63Jp0o+EZG8KUGVQK32hRIR6TMlqBKo1QhKRKTPlKBKQDvrioj0nRJUCdQ0jyCI15Fs3Upyx/aowxERqQhlv+U7gJn9GHgP0JnRfIS7v5Kl7ynAt4DJwALgYndfUpJAcwiCgPjIvelY+xodm1ZTP25SlOGIiFSEikhQwDTgXHd/qKdOZrYX8BvgUuAe4BPA78xsirsnix9mbrUjlKBERPqi7Kf4zKwROASYn0f384AX3P3X7t7h7jcA9cBpxYwxH1rVXESkb8piBGVmdcCoLE+lgEmEU3vfN7PjgNeAz7v7fVn6HwYs6tbmwOHAw4WLuO+0qrmISN+UywhqOrAyy+N1YCjwJPBFYAIwG/ilmU3Ncp5moKVbWwvQVJyw86dVzUVE+qYsRlDu/jgQ9NBlbsb3vzazfwPOJiyCyLSdPZNRE7BtoDEOlNbjExHpm3IZQeVkZmeZ2Qe6NdcBO7J0XwRYt7ZD2HPar+S6tt3o3LxW226IiOShLEZQvagBvmVmLwLPAO8lnBL8YJa+dwNfN7P3pL+/EkgCj5cm1Nxi8Tpqho4isXUDnVvW7ZzyExGR7PJKUGZ2OHAGcDQwFkgAq4B5wH3uvrhYAbr7b8zsGuDnwDjgJeBd7v5qOrZbgAPc/Qx3X2NmZxHeB/VD4AXgLHdvL1Z8fVE7Yu8wQW1crQQlItKLHhOUmZ1IWJwwHXiKcKpsMeGoZi/g/YQjlieA69z9D8UI0t2/C3w3x3OXdzv+A3BkMeIYqPjIcfDai3RsXEXjxCOiDkdEpKzlTFBm9n/AG4CbgRnuvilHv2HARcBNZvacu19SjEAHA61qLiKSv55GUPe7+6W9ncDdtwC3Area2QUFi2wQ6qrk0826IiK9y1nF5+6/7uvJ3P1XAwtncKsdOR6Ajg0rI45ERKT85Vsk0QT8P+An7v6ymd0KzCS8LjXT3fUbNw+1o9IJauNKUqkkQVD2Vf4iIpHJ9zfktwgLIurM7BzgA8BVQCvw7SLFNujUNDYTaxpGqqONxNaNUYcjIlLW8k1Q5wAXufsLwLuBh939+8B/Am8vVnCDUe2oCQB0bFgRcSQiIuUt3wTVCKw2sxjwL0DXthcpwnuiJE87E9R6JSgRkZ7ku5LEPMJrUGuBkcDdZjYBuA74a5FiG5TqRmsEJSKSj3xHUB8jvFn3SuDf3X0F8BnCde8+XqTYBqWuEVS7RlAiIj3q6Ubd6cBf3T3p7ouAN3Xr8ll331rU6Aah2tFdpeZKUCIiPelpiu9GwMzsScLtLuZmrrmn5NQ/8ZHjgIDOTWtIJToIamqjDklEpCzlTFDuPj29jNFpwOnAJ82shnBn2rnAo+6uWuk+isXriA8fQ+fmNXRsWkPd6H2iDklEpCz1WCSRXsbo7vQDM5tIWFZ+IXCLmS0hLDn/XLEDHUxqR48PE9T6FUpQIiI59Gk/KHf/J7vW3YsBxxCOrqQPakdNoPWVBVrySESkB3knKDM7mXB18/puT7UUMqBqoJt1RUR6l+9afDcB/wG8yp5braeAbxY4rkFt55p8SlAiIjnlO4L6V+BSd7+9mMFUi9qum3XXa4pPRCSXfG/UbSFcuVwKID5sL6iJk9i2gWR7a9ThiIiUpXwT1JeBG9NVfDJAQayG2pHjAO0NJSKSS75TfC8CXwEWm9keT7p7TSGDymRmtxBu9ZFpCHCNu38lS/8vAbOA9ozmd7n748WKsT9qR02gY91yOjaspH7cpKjDEREpO/kmqO8RLgp7GyWu2nP3y4HLu47N7ArCtQFvzvGSacDH3f2WEoTXb7uuQ70ecSQiIuUp3wS1H3CGu79SzGB6Y2aTga8DJ6RvIs7mKMIpybK2q5JPU3wiItnkm6AeBk4EipKgzKwOGJXlqZS7r844/gbwfXdfkOM844FxwCwzOx5YD9zg7rcVOuaB6lpBQvtCiYhkl2+C+htws5mdDywGOjKfdPdPDzCO6cBjWdoTpGM0szcSrlrR0wWbscAThNN/706f9x4zW+Pu9w8wxoLKvFk3lUoRBEHEEYmIlJd8E9TphJsWNrPnthupgQaRLmDo7Tf0pcBvuo2oup9nAXByRtMTZnYHcB5QVgkq1jSMWH0TybYWki1bqBkyPOqQRETKSl4Jyt1PKXYgeTgH+ERPHczsBOBod78po7mOPVe/iFwQBNSOmkDbysV0bFihBCUi0k3O+6DM7Atm1pjvicxsaLrEu+DMbAzh1N6feunaClxvZmeaWczM3ga8DyjLFTC6Kvm0u66IyJ56ulF3M/CCmX3dzI7L1sHMAjM7xsy+BSwCNhUjSOBAoM3dN2SJ4RYzexDA3Z8BLias9NsKfAe4xN3LchUMrcknIpJbTxsW3mRmvwI+Dcw1s07CG3bXEV4vGkO4unkA/Ah4i7u/Wowg3X0e0JDjucu7Hd8J3FmMOAptV6GESs1FRLrrbcPC5cDHzewzhMUHRwF7A0nCyr4vAo+5e1uR4xyUtO2GiEhu+RZJbCesgiurSrhKVzt61826qUQnQU2f9o8UERnU8l0sVoogVtdIfPhYSHTSsXFV1OGIiJQVJaiI1Y3ZD4D2tUW5fCciUrGUoCJWN3Z/ANrXKEGJiGTKK0GZ2YFFjqNq1Y5JJyiNoEREdpPvVfnFZvZX4CfAne6+vogxVZW6dILqUIISEdlNvlN8k4D7gCuAFWZ2r5ld2JeVJiS7utH7QBCjY8Mqkh2q1hcR6ZJXgnL3V939enefSngv1ALgM8BqM7s9vaSQ9EMQr00veZSiY502LxQR6dKfIonlwBLCvaHihKub/9jMPL0Hk/SRKvlERPaUb5HEEDN7n5ndA6wCrgVeAo5Jj6r2JdzP6RfFCnQwq1OhhIjIHvItklgDtAN3Ae9I79+0k7snzWwu8NbChlcdlKBERPaUb4K6BLgn25p7ZjbW3de4+12ECUz6SPdCiYjsKd9rUHOAYd0bzWx/wmtRMgDxEXsTxOtIbF1PYsf2qMMRESkLOUdQZnYRMCN9GAA/MLPuI6gDgD32aJK+CWI11O61L+2rXqFj7WvU7HdI1CGJiESupxHUw8A2oOtP+tb0912PbYRbbpxbzACrha5DiYjsrqcNC9cBlwKY2VLgBndvKU1Y1Uel5iIiu+tpiu+dwMPu3gHMA042s6x93f2B4oRXPerGHgCoUEJEpEtPVXz3AeMIS8zv66FfCqgpZFDVKHOKL5VKEQRBxBGJiESrpym+WLbvi83MrgJOcvdzM9r2B34IHEeYMP8j16jNzEYAPwBOJ7xO9jl3v63ogQ9QzdBRxOqbSLZuJbF9E/HmkVGHJCISqbwTj5ldamYXZBz/0swuLlQgZtZsZjcA38jy9BzgOWA08CFgjplNynGq/wUSwHjgTOB6MzupUHEWSxAE2npDRCRDvksdXQPcyO5Tec8DN5nZJwoUy/3ARODWbu89BTga+Ly7t7v774F7gMuyxNkEXAD8l7u3uPt84PvAhwsUY1Ht2nrjtYgjERGJXr4rSXwEuNDd53Y1uPt1ZvY08D/ATb2dwMzqgFFZnkq5+2rgIndfYWbXEo5+uhwGvOrumXewvgQcm+VcUwivib3cre/ZvcVXDrSihIjILvkmqJHAsiztS4C98zzHdMIFZbtLAHF3X5Hjdc1A9/L2FqApR98d7p7Ko2/Z0b1QIiK75Jug/grMMrMPuXsngJnVAFcTlqD3Kr3AbH9K07YD3TdGbCIsgMjWt8HMgowklatv2dl1L9RrpFJJgqBktSkiImUn3wT1KeBR4FUze45wGu3w9OvPKFJsXRYB+5tZo7u3ptsOSbd39zJhEpzIrjUCc/UtOzVNw6gZMoLE9k10blpD7chxUYckIhKZfHfUXQAYMBtYDLwIfBk42N2fLV544O5OuIPvbDOrN7NTgHOAn2Xpuw24G/hquipwKmHV3x3FjLGQ6vaeCEDbKq3BKyLVLe85JHdfD/wOmAv8AXjc3bcWK7BuzgcOJbwH6gfAZe6+EMDM3mpm29L3SkFY0JEkvGb2ADDb3R8sUZwDVj9hMgBtKxZHHImISLTymuIzs2bCG2UvADoIp9HiZvYwcH63CrsBcfdrs7S9Ro6pRHd/krA4out4I3BRoeIptfoJBwPQtlIJSkSqW74jqG8SXnM6nrBgoSH9/QTga8UJrTrVjz8IgLaVS0glExFHIyISnXyLJM4DZrj7UxltT5nZR4FfAR8reGRVKt48gpphe5HYso6O9St2VvaJiFSbfEdQMWBdlvYNZEyvSWE0TOgaRWmaT0SqV74J6g/AtenVIAAws3rgC8CTxQismu2c5lOhhIhUsb7cB/VH4DUzm59umwrsAN5RjMCqWf0EJSgRkXzvg1pMWOZ9HeHyRguBzwGHuvtLxQuvOtWPmwQEtK1eSqqzI+pwREQike8Iqqt8++YixiJpsYYh1I6eQMf612lbs2znNSkRkWrS05bv8wiXNOqVu2dbWVwGoH7CQWGCWrFYCUpEqlJvW75LROrHH8S255+gbeXL6DKfiFSjnrZ8/2IpA5HdqVBCRKpd3tegzOw9wH8CBwPTgH8HVrn7jUWKrarV7X0gxGroWPc6ybZWYvXddxwRERnc8t3y/RLCnXPvArruhXoJ+LyZzSpOaNUtFq+jbuyBQIq2VUuiDkdEpOTyvVH3auAKd/8q4Q64uPsPgH8jXD1cikArm4tINcs3QU0Gns7SPh/QrnpFsmvhWCUoEak++SYoB96Wpf09hFN9UgQNXVtvaAQlIlUo3yKJzwK/MrOj06+53MwOAt5FuEeUFEHtXvsQ1DbQuXktie2bqRkyPOqQRERKJt+ljh4EjgXqCZc5Op1wHb7j3P2e4oVX3YJYDfXj0lvAa5pPRKpMTytJvBN4yN2TAO7+AnBJieKStPp9DmbHay+y47WXaDroqKjDEREpmZ6m+H4LrDeznwO3u/v8HvoWjJldBZzk7udmtB0F/DdwBLAF+AFwnbvvsRSTmU0CFgMtGc1z3P2DRQ28SBoPeCOb/3oPrUufjzoUEZGS6ilB7QNcmH5caWYvALcDP3P3FYUOxMyaCfeXuhq4J6O9Cbgf+DJwCjAJ+B2wCvhellNNA55y9+MKHWMUGvY7DGI1tK1cQnLHdmINQ6IOSUSkJHJeg3L3Ne7+bXefTpgUfgZcDCwzs7lm9n4zK+TyBvcDE4Fbu7XvB/zF3W9294S7vwz8Bjghx3mOIix/HxRi9Y3hskepJK2vLoo6HBGRksmris/dlwJfBb5qZocB7yOs7PsfM/u1u/9bb+dI78Y7KstTKXdfDVzk7ivM7FpgfMZ7OzCj23nOIPvoCcIR1BAz+wfhdvQPAJ9y9029/qBlqvGAw2lb7rQuW8iQKcdEHY6ISEnkex/UTu6+CLiRMGG9TDiqysd0YGWWx+vp8/Y6bZjeZv7nhNeXbsnRbSPwMHAMYbLan9zJrCI0TjwcgB26DiUiVaQvi8UOB84lvDn3NOAV4KfAefm83t0fB4K+h7jz/ccBvwaSwNvcvTXH+1yYcbjZzD4L/NHM4u7e2d/3j1L9PlMI4nW0r1mm+6FEpGr0mKAyktK7CVeS2ALMAa5193nFD29nHIcRFkb8Hviwu7fl6NcEXAt8Iz1tCOHitp2k1xCsRLF4HQ37HULrP5+jddlCmg97S9QhiYgUXU/3Qd1POFJKAvcSrhjxUKlHIWY2EphLWCr+qZ76unuLmZ0OjDazjwEjgOuBH2UrSa8kDQccHiaopc8rQYlIVehpBNUMfBS40923lCiebC4mLHm/wswuz2i/190vMrO3Ag8Ch7n7q4QFFd8BVhCOmuYAPSa2StB44OFsBN0PJSJVo6cddU8qZSAZ73ttt+NvA9/uof+ThMm063gpcFaRwotM/fhJBPVNdG5cRefmtcSHj4k6JBGRoupzFZ9EI4jV0Lj/YYBGUSJSHZSgKkjjgWG5eeuyhRFHIiJSfEpQFWRnglr6PKlURdd8iIj0SgmqgtSO2Z9Y0zASWzfQsaHgyyGKiJQVJagKEgTBzlGUVpUQkcFOCarCdCWollcGzXq4IiJZKUFVmKbJ0wBofWUByfasqz2JiAwKSlAVJj5sNPX7HkKqs52Wl5+JOhwRkaJRgqpAzYceD8C2F/8ccSQiIsWjBFWBhhwSJqjWJc9qmk9EBi0lqAoUHzaa+n0snOZb/PeowxERKQolqAo1JD3Nt/3Fv0QciYhIcShBVajmQ44DoGXxMyTbd0QcjYhI4SlBVaj48DHUTzg4nOZb8mzU4YiIFJwSVAUbcuh0ALarmk9EBiElqAo25NCMab6OtoijEREpLCWoClY7fGw4zdfRRssSVfOJyOCiBFXhhqSLJVTNJyKDjRJUhesqN2/5xzwSrdsijkZEpHDiUQfQnZldBZzk7udmtJ0KPAxkLpvwNXe/LsvrRwA/AE4HtgGfc/fbiht1dGpH7E3jxKm0/nMBWxc8yojjzok6JBGRgiibBGVmzcAXgKuBe7o9PQ24090vzONU/wskgPHAFOB3ZvaKuz9RyHjLyfBjzqT1nwvY8vSDDD/2XQSxmqhDEhEZsHKa4rsfmAjcmuW5o4BeN0AysybgAuC/3L3F3ecD3wc+XMhAy03jQUcSHzmOzs1raXn56ajDEREpiJKNoMysDhiV5amUu68GLnL3FWZ2LeHoJ9M0YIyZXQEEwC8Ip+6611ZPAVLAyxltLwFnF+BHKFtBEGP40Wew/uHb2DzvAYbYm6MOSURkwEo5gpoOrMzyeB3A3Vdke5GZxYHlwN3AocCpwNuAPa4/Ac3ADndPZbS1AE2F+RHK19AjTiGoa2DHsoW0rV4adTgiIgNWshGUuz9OOPrp6+s6gdMymhab2Wzga8Cnu3XfDjSYWZCRpJoIiyUGtVjDEIYecQpbnn6QLU8/yJgzr4g6JBGRASmna1BZmdk+ZnZjeoqwSx2QbYXUlwmT4MSMtkOARUUMsWwMO/qdAGxb+AcSLVsjjkZEZGDKPkEB64GZwOfMLG5mBwOfA/YoHXf3bYRTgV81s2Yzmwp8CLijlAFHpW70BBonH0mqs52t8x+JOhwRkQEp+wTl7juAM4ATCZPVH4A7gW8CmNlbzWybme2ffslHgCSwDHgAmO3uD5Y88IgMP+ZMADY/8xCpRGfE0YiI9F/Z3AfVxd2vzdI2Hzg5R/8nCYsjuo43AhcVKbyy1zhpKrWjJ9CxfgVbnn2E4Ue/I+qQRET6pexHUNI3QRBj1MkzAdj4hzkkdmyPOCIRkf5RghqEmuzNNOz/BpKtW9n0xzujDkdEpF+UoAahIAgY/bZLgIDN8x6kY0PWW8xERMqaEtQgVT9+Es1HnALJTtY/+uOowxER6TMlqEFs1MnvI6htoOUf82hd+nzU4YiI9IkS1CAWHzqSEdNnALD+4dtIJRMRRyQikj8lqEFu+JvPIj5sL9rXLGPzU/dHHY6ISN6UoAa5WG09o//lgwBseOwn7Fj+UsQRiYjkRwmqCgyZcgzD33wWJBOsvuubJFq2RB2SiEivlKCqxKhT3k/9PkZi63rW/PbbpFLJqEMSEemRElSVCGri7H3eJ4k1DqX1lWfZ9Oe7ow5JRKRHSlBVJD5sL8aecyUAG5+YQ+s/n4s4IhGR3JSgqkzT5CMZ8ZbzIZVk1Z1fo3XZwqhDEhHJSgmqCo088b00H34yqY4drJozm5Z/Log6JBGRPShBVaEgVsOYd/07Q6eeRqqzndW/vJ6WJc9GHZaIyG6UoKpUEKthrzMvZ+i0t5PqbGfVndez3Z+KOiwRkZ2UoKpYEMTY6x0fZtjRZ0Cik9W/+hrrH/0xqURH1KGJiChBVbsgCBj99ssYedJFEMTY/NffsuL2a+jYsDLq0ESkyilBCUEQMPKEC5jwr9cRHz6GtpVLWP7DT7F1we91Q6+IRCYedQDdmdlVwEnufm76+K3Ag9261QP/dPcpWV4/CVgMtGQ0z3H3DxYp5EGjYd9D2OeD32Ddg7eyfdGfWHvfd9n89EOMOmUmTZOmRh2eiFSZsklQZtYMfAG4Grinq93dnwSaM/rtB8wDPpbjVNOAp9z9uOJFO3jVNAxh7LlXsW3yNDY89hPaVy1h1c+/ROOBhzPy5Jk07HNw1CGKSJUomwQF3A+sBW4FxvfQ7/+AO9x9bo7njwLmFzi2qhIEAUOPOJkhhx7PlnkPsOkvd9O69HlafzSL+n2NYW96G0MOPZ5YXWPUoYrIIFayBGVmdcCoLE+l3H01cJG7rzCza8mRoMzsXOAw4Owe3moaMMTM/kE48noA+JS7bxpI/NUoVlvPiOkzGHrk6Wz6y91seeZ3tC131i531s39Ic2HncCQQ6fTuP9hBPHaqMMVkUGmlCOo6cBjWdoTQNzdV+RxjmuA6929tYc+G4E/AzcBjcCPge8B7+lbuNKlprGZ0adezMgT3s32F//MlvmP0rb8JbbOf4St8x8hqK2n8YA30jh5Go0HvpHa0RMIAtXfiMjAlCxBufvjQNDf15vZEcAbgNt7eZ8LMw43m9lngT+aWdzdO/v7/gKxugaGTj2VoVNPpX3dcrY9/wQti/9O+8FInp0AAAk9SURBVJqltCx+hpbFzwAQ1DdRP24S9RMOon7cJGpH70PtqPHEausj/glEpJKU0zWo3pwDPOjuOXfbM7Mm4FrgG+lpQ4A6oJNwpCYFUrfXvow6ZSajTplJ59YNtCx5ltZXnmXH8n+Q2LqeHcsWsqPbQrQ1w/aibtR4aoaNIT5sFPGho4kPG01N03BiQ4aFX5XERCStkhLUccCjPXVw9xYzOx0YbWYfA0YA1wM/cvdUCWKsSvGhoxj2ptMY9qbTAOjcupG2lYtpW7GY9rXL6Fi/go6Nq0hsWUfrlnU9niuI1xFrbCZW35TxaCSobSBWW09QW0+stoEgXpt+1IVfa+JQEyeIxQlq4gSxGojVENTUQFBDEItBENv5lVgsnIYMgl3fE4TH6a9BkB7w72xj13N7tJPx+mw/2K72oPtEwh4v6cdEQ673FSmBIFZTlPNWUoI6ENjjOlXGfVKHufurwAzgO+m+CWAO8KnShSnxoSOJDz2GIVOO2dmWSibo3LSajg2r6Ny6ns4t4SOxdT2Jli0kWjaTaNlCqrOdxNYNJLZuiPAnEJG8BTFGnTKTEcefW/BTl12Ccvdrc7S/IUf7bvdJuftS4KxixCb9F8RqqB01gdpRE3L2SaVSpNp3kNyxjWRbS/jY0UKyvZVkRxupjjZSHTtItreRSnSQ6uwg1dkePpIJUolOSHSSSnaSSibT3ychmSCVTEAqGa6MkUyGxwDJ5K72FJBeOWPnChqpFJAKnyMFqVT628z2Xc/l+MF2fUu3Pnu8pB8D/VzvK1IKsVg4g1EEZZegpHoFQUBQ30isXvdXiYjW4hMRkTKlBCUiImVJCUpERMqSEpSIiJQlJSgRESlLSlAiIlKWlKBERKQs6T6o3dUArFq1Kuo4RESqQsbv2z3WS1KC2t14gJkzZ0Ydh4hItRkPLMlsUILa3TzgrcBKtPq5iEgp1BAmp3ndnwhSWsdLRETKkIokRESkLClBiYhIWVKCEhGRsqQEJSIiZUkJSkREypISlIiIlCUlKBERKUtKUCIiUpa0kkSBmNlU4BbgCOAV4FJ33+PO6GpkZpcCtwJtGc0fdffbIwopcmZ2LHCfu49NH9cBNwMXEK5i8k13/2qEIUYmy79NPbAVaM/o9md3f3sU8ZWamZ0OXA8cDKwBbnD3W6vhM6MEVQDpD8pvgZuAE4HzgblmdoC7b4k0uPIwDfiGu8+KOpComVkAXAbc2O2pLwIGTAaGAw+Z2evu/uMShxiZHv5tDgc2uPu40kcVLTPbD/g18AHC3zFHAb8zs6XAyQzyz4ym+ArjZKDW3W9y9w53nwO8ALw32rDKxlHA/KiDKBNfBK4Avtyt/QPAbHff6O5LCX9Jf6TEsUUt179NNX9+DgR+5u53u3syPSvzOPAWquAzowRVGIcBL3Zre4nwL7+qZmY1hNOeF5vZCjNbbGaz0n8tV6Nb3P0o4OmuBjMbQbhY5qKMftX4+dnj3yZtGjDWzJ4zs9VmdqeZ7RNBfCXn7k+6++Vdx2Y2inBB62epgs+MElRhNAMt3dpagKYIYik3Ywh/4dwOTCScL78i/ag67r4iS3Nz+mvmZ6jqPj85/m0AtgN/Ak4jnNJqBe4uVVzlwsyGA/cAfwOeSTcP6s+MrkEVxnagsVtbE7AtgljKiruvAk7KaJpvZt8hvE73P9FEVXa2p79mfob0+Ulz909mHpvZJ4G1Zrafu78WUVglZWZTCK9BLQJmsuuzMqg/MxpBFcYiwr/sMh3C7sPvqmRmbzCzL3ZrrgN2RBFPOXL3jcAqdv8M6fOTZmZfMrNDM5rq0l+r4jNkZicSjpp+A1zg7juq5TOjEVRhPAYEZnYVYdnn+YTXXapuGiKLTcDVZrYc+CFwJPBx4GORRlV+7gC+YGbPEU75fQr4VrQhlY0jgKPN7H3p428B97v72ghjKgkzmwzcB1zj7t/p9vSg/8xoBFUA7t4OnEGYmDYA1wDnVsN/oN64++vA2YTVRVsIS2avc/dfRRpY+fk8sJCw+nMe4b/TLZFGVD4uAzYCi4GlhPdDXRxlQCX0UWAo8FUz25bx+BpV8JnRjroiIlKWNIISEZGypAQlIiJlSQlKRETKkhKUiIiUJSUoEREpS0pQIiJSlnSjrkiRmdmPCFeezuWLhCtUPwYMdfeSLFeTXsj3T8C/uvs/eugXA/4KXOzuXorYREAjKJFSuJJw5enxhFuzAByb0XYj8Of099uzvL5YPg4s6Ck5Abh7EvgSg+wmUCl/ulFXpITM7I3A88DE9B4+UcXRALwKnOruC/N8zRLgMnd/vJixiXTRFJ9IGTCzk8mY4jOzFHAR8BnCBUGfBt4P/CfhMj9bgM+4+x3p1w8FvkG4nUkK+D1wZQ9bWFwIbMpMTmb2X8CHCbdIeRH4rLs/mPGauwlHg48X4EcW6ZWm+ETK1/XAJ4DjgP2BvxMmpmOAu4BbzaxrL6nvESayfyHc3iRFuDV4rj9CzwQe6jowsxnp93o/4arY9wN3mtmwjNc8BLyth3OKFJQSlEj5+q67P+bu8wlXtN5GOKpx4JuEewFNNLNJhCOi97n7vPSo6GLC7cLfkePcRxMuMtrlQKANWJaeevwScB7QkdFnEeGq2YcU5KcT6YX+EhIpX4szvm8Blrp710Xjrr2Q6oED0t+72W7bkjURjqruy3LuvYF1Gcc/Iaw0fMXMniHcufU2d2/N6LM+/XVsH38OkX7RCEqkfHV0O07m6BdP9z0SeFPGYwpwW47XJIGg6yC9NcxRhCOuPwOXAM+lizq6dP2+SOT9E4gMgBKUSOV7EagFhrj7YndfDKwEbiBMUtmsIiyGAMDMzgM+4u5z3f1KwpHXVuCdGa8Zk/FakaLTFJ9IhXN3N7N7gB+b2UeBtcBswuKKl3K87BlgasZxDXCDma0mrBg8DhiX/r7LVHZtHChSdBpBiQwOHyBMJr8h3F11OHC6u2/K0f9+wmo/ANz9TuALhKOufwBfBj7m7r/PeM2JwEPurik+KQndqCtShcysiXD79He4+9/z6B8DlhFWCj5Z5PBEAI2gRKqSu7cQjpY+mudLzgFeUXKSUlKCEqle/w0cYd1q07tLj56uAS4vSVQiaZriExGRsqQRlIiIlCUlKBERKUtKUCIiUpaUoEREpCwpQYmISFn6/8XmJ6wJ+7iTAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def plot_velocity(results):\n",
    "    plot(results.v, color='C1', label='v')\n",
    "        \n",
    "    decorate(xlabel='Time (s)',\n",
    "             ylabel='Velocity (m/s)')\n",
    "    \n",
    "plot_velocity(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From an initial velocity of 0, the penny accelerates downward until it reaches terminal velocity; after that, velocity is constant."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run the simulation with an initial velocity, downward, that exceeds the penny's terminal velocity.  Hint: You can create a new `Params` object based on an existing one, like this:\n",
    "\n",
    "`params2 = Params(params, v_init=-30 * m/s)`\n",
    "\n",
    "What do you expect to happen?  Plot velocity and position as a function of time, and see if they are consistent with your prediction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>height</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>-30.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.0025 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.019 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>18.0 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "height                      381 meter\n",
       "v_init           -30.0 meter / second\n",
       "g             9.8 meter / second ** 2\n",
       "mass                  0.0025 kilogram\n",
       "diameter                  0.019 meter\n",
       "rho         1.2 kilogram / meter ** 3\n",
       "v_term            18.0 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "v_init = -30 * m / s\n",
    "params2 = Params(params, v_init=v_init)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'A termination event occurred.'"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "system2 = make_system(params2)\n",
    "results, details = run_ode_solver(system2, slope_func, events=event_func, max_step=0.5)\n",
    "details.message"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3dd3xUxfrH8c9mk5CE3hOaIGXoqHRQWixg73oVBUWvIHaqeq+KPxQQBEGavWC9itiwoKGJ9CYlMPQmhCY91CS/P85GlwhxgWTPbvJ9v155ZXfmnN0nuObJPGfOjCcjIwMREZFQE+F2ACIiIqeiBCUiIiFJCUpEREKSEpSIiIQkJSgREQlJSlAiIhKSIoP9hsaYYsAS4Blr7bvGmGhgJHAzkAYMtdYO8Dv+VuBFIAGYBnS21u4IdtwiIhJcboygxgLl/Z73AwxQFWgMdDLG3A1gjKkNvAV0BkoCq4FPghmsiIi4I6gJyhjTCSgCLPVr7gS8YK3dY63dAAwBHvD1dQS+sdbOsNYeAZ4EWhpjqgcxbBERcUHQSnzGmCrAs0AL4AdfWzGc0l2y36ErgXq+x7WB+Zkd1tpUY8xmX//qAN+3AM7IbBtOCVFEREKHFycPzLPWHvXvCEqCMsZ4gQ+AntbaFGNMZlch3/dUv8NTgTi/fv++rP2BaAz8ckYBi4hIsF0CzPBvCNYI6r+AtdZ+kaX9kO97rF9bHHDQrz+Wk/n3B2IbwIcffkh8fPwZnCYiIrktJSWFO++8E3y/q/0FK0HdDpQzxtzoe14YGA00AVJwJkn87uuryV8lv2RfHwDGmDigEieXBP9JGkB8fDwVKlQ42/hFRCR3/e0STFASlLW2pv9zY8xi4BXfNPODwLPGmCU4Jb2ewHDfoR8BM4wxbYBZwABgkbV2VTDiFhER94TCjbrPAMuA5cA8YDzOVHSstUuBe33PdwF1gFvcCVNERIIp6DfqAlhrL/B7fATo7vs61bHjcZKWiIjkI6EwghIREfkbJSgREQlJSlAByMjIcDsEEZF8RwnqH3zwwwru7f8Ty9bucjsUEZF8RQnqH5w4kc6uvYd59vVZzFn2t/vIREQklyhB/YO7rqxN++aVOXYinRffm8fPcze6HZKISL6gBPUPvBEeHrypPrdfZkhPz2D4p4sZPzmgdWpFROQcKEEFwOPxcGf7mvz7emeR9XcnJvPW18tIS9fkCRGR3KIEdQauueR8et7ZEG+Ehy+nraX/23M4dPi422GJiORJSlBnqPVFFXj+geYUjoti/ort9Bwxna07z2RxdRERCYQS1FmoX600Qx9rzXnxhdmy4yBPDJ/OQrvD7bBERPIUJaizFF+yIC89fAnN6sZz6PBx+r0xiy+mrNZNvSIiOUQJ6hzExUTxZKcm3HZZDdIz4J1vk+n/9lwOpB5zOzQRkbCnBHWOIiI8dGxfi//c04SCsVHMTU7h0aFTsRv/cDs0EZGwpgSVQ5rWTWD4E22oUakYO/ccpu+oGXw1fa1KfiIiZ0kJKgeVLRHHwO6XcG2r8zmRlsGbXy3j+bfmsPfAUbdDExEJO0pQOSwqMoL7r6vHU50bUzDWmYr+8JApzF+x3e3QRETCihJULmlerxyv9mhLvaql2HvwKP3enM3YL5Zw9Hia26GJiIQFJahcVLp4LP27tuCeq2sT6fUw8df1PD5sKms273U7NBGRkKcElcsiIjzc2LY6Qx5pRYUyhdi8/SA9Rkznwx9WcvxEutvhiYiELCWoIKlaoRjDHm/Nta3OJz09g09+svQcPp0N2/a7HZqISEhSggqimOhI7r+uHi8+2JKyJeJYt3Ufjw+bymdJq0hL02hKRMSfEpQL6lUtxYgebWjfvDIn0jJ4/7sV9BwxnfVb97kdmohIyFCCcklcTBTdb27A8/9uTunisazZso/Hh03jox91bUpEBJSgXHehKcPInm25skVl0tIz+HiS5YlXprFq0x63QxMRcZUSVAiIi4mi200NePHBliSULMiGbfvpNWI6b329jCNHT7gdnoiIK5SgQki9qqUY0bMN17euCsCX09bSfcgU7TUlIvmSElSIiYmOpMu1dRnyaCuqlCvCjj9Sefb1WQz7eCH7D2kbDxHJP5SgQlT1isUZ+lhr7r6yFlGREUyev5lug5KYPH+zVkgXkXxBCSqERXojuCWxBiN7tqV+tVLsP3SMYR8v5JnXZrF110G3wxMRyVVKUGGgXOlC9O/agkdvu5DCcVEsXr2ThwdP4X8/r9KUdBHJs5SgwoTH4+HSJpUY0yeRtg0rcOxEOuO+X8Fjw6ayfN1ut8MTEclxSlBhpmihAjxxR0P6P9CChFIF2ZRygL6jZjDi00UcSNUkChHJO5SgwlSDGqUZ2bMtt11Wg0ivh5/mbqLrwCQmz9+kSRQikicoQYWx6CgvHdvXYkSPttStWtI3iWIRT4+ZyebtB9wOT0TknChB5QEVyxbmxW4tefS2CylSMJqla3fxyMtT+OD7FdrBV0TClhJUHuE/ieLypudxIi2DT39exUODJ7Ng5Xa3wxMROWNKUHlMkYLRPHzrBQx66GIqJxQhZXcqz70xm4HvzWPX3sNuhyciEjAlqDyqdpWSDHu8NfdeU4eYaC+/LtnKgy8l8eW0NZzQ5ogiEgaUoPKwSG8EN7SpxujeiTSvl8Dho2m89fVyHh82jeT1undKREKbElQ+ULp4LE91bsKz9zWjbIk4NmzbT5+RMxj+ySL2HTzqdngiIqekBJWPNKpVllG923HbpTWI9Ebw8zzn3qkfZm0gPV33TolIaFGCymcKRHnp2KEWI3u15YLqpTl4+DijPv+NXq9OZ82WvW6HJyLyp8hgvpkx5mrgRaAKsAN4yVr7mjEmGhgJ3AykAUOttQP8zrvVd14CMA3obK128TsX5UsX4vkHmjPjt628+dVSVm3aS49XptGhRRU6dqhFodgot0MUkXwuaCMoY0wC8DnQx1pbGLgFeMUYcxHQDzBAVaAx0MkYc7fvvNrAW0BnoCSwGvgkWHHnZR6Ph0suKM+YPolc16oqeDxM/HU93QZq3ykRcV/QEpS1dhtQ2lr7vTEmAifZnAAOAJ2AF6y1e6y1G4AhwAO+UzsC31hrZ1hrjwBPAi2NMdWDFXteFxcTxX3X1eWVx1tTu0oJ9h48yrCPF/Lk6F/ZuG2/2+GJSD4V1GtQ1toDxpg44CgwCRgF7MQp3SX7HboSqOd7XNu/z1qbCmz265ccUqVcUQZ2v5jHbr+QooWiWb5uN48OncpbXy8j9chxt8MTkXwmqNegfI4ABYH6wHdA5vIGqX7HpAJxvseFsvRl7Zcc5PF4SGxciaZ14hn3/Qq+n7WBL6et5ZfFv9Pl2rpc3KAcHo/H7TBFJB8I+iw+a226tfaYtXY+8DrQyNcV63dYHJC5p/mhLH1Z+yUXFIqLpttNDRj6aGtqVCrG7n1HeGncfJ55fRZbdmildBHJfcGcJNHaGLMgS3MBYA+QgjNJIlNN/irrJfv3+UqElTi5JCi5pFrFYgx+uBXdb25AodgoFq/aycNDpvD+d8kcOXbC7fBEJA8LZolvMVDeGPMEMBxoCnQBbsBJUM8aY5bglPR6+o4B+AiYYYxpA8wCBgCLrLWrghh7vhYR4aF988o0r5fAexOT+WnuJj5LWs3UhVu4/7p6NKsbr7KfiOS4YM7i2wdcCdwI/IFT3rvPWjsNeAZYBiwH5gHjgbG+85YC9/qe7wLq4ExRlyArWqgAj9x2IYMfvoTzyxVl557DvPjuXPq9OZttuw65HZ6I5DGevH6vizGmMrA+KSmJChUquB1OnpGWls73szYw7vsVpB45QVRkBDe3q85N7apTIMrrdngiEia2bNlCYmIiQBXfbUZ/0lJHcla83giuvvh8xvZNpG3DChw/kc7HkywPDZ7MvOQUt8MTkTxACUrOSfHCMTxxR0MGPNiS8+ILk7I7leffmkP/t+ew/Y+sdweIiAROCUpyRN2qpXjliTZ0ubYusQUimbM8hQcHJfHpT5Zjx9PcDk9EwpASlOSYSG8E17euypg+7Wh9YQWOnUjngx9W8tCQKSxYud3t8EQkzChBSY4rWTSWnh0b8mK3llQsW5htuw7x3BuzefHduexQ2U9EAqQEJbmmXrVSjOjRhnuurk1MtJdZS7fR7aXJ/O/nVRw/obKfiGRPCUpyVaQ3ghvbVmdMn0QuuaA8x46nMe77FTw0eAoLtaWXiGRDCUqColSxWHrf1Yj+D7SgQplCbN11iGdfn8XA9+axc8/hf34BEcl3lKAkqBrUKM2IHm25+8paFIj28uuSrXR7KYnPJ6/m+Il0t8MTkRCiBCVBFxUZwS2JNRjTO5GW9ctx9Fga701M5pGXp/Dbqp1uhyciIUIJSlxTungsfTs1pt+/m1O+dEG27DjIf16byaD357F7n8p+IvmdEpS47iJThld7OmW/6CgvM37bSteBSXwxRWU/kfxMCUpCQlSk1yn79WlH83oJHDmWxjvfJvPo0Cn8tlplP5H8SAlKQkqZ4nE81bkJz93fjIRSBdm8/SD/GTuTwePmq+wnks8oQUlIalizLKN6taVj+5pER3mZvvh3ug1KYsLUNZxIU9lPJD8IeEddY0xZoCFQBkjD2QV3obV2dy7FJvlcVKSX2y4ztGlYkTe/WsrsZSm8/c1yfp63ia431KdetVJuhygiuSjbBGWMiQTuAB4DGgDHgD2AFyjhO2YOMBr4xFqrP20lx5UtEcfT9zRl/ortvD5hKZtSDvDUmF9pfWEF7rmmNiWLxrodoojkgtOW+IwxrYElwN3AW0ANIM5aW85aWxaIBi4EPgIeAlYaY9rkesSSbzWqVZaRvdpyZ/uaREdGMG3RFroNmsyX09aq7CeSB2U3guoB3GatXXqqTmttBrDM9zXaGHMh8DwwNaeDFMkUHeXl9ssMbS6qwJtfLWPO8hTe+noZP8/dSNcb61O3qsp+InmFJyMjw+0YcpUxpjKwPikpiQoVKrgdjuSweckpvP7lUlJ2O9t4tGlYgXuurkOJIjEuRyYigdiyZQuJiYkAVay1G/z7zmSSRBxQBSiQtc9au/AcYxQ5K41rx9OgemnGT1nD50mrmLpgC3OWpXBn+5pc3bIKXq8mqoqEq4ASlDGmIzAWiAU8WbozcCZNiLgiOsrLvy43tG1YgTe+XMbc5BTe/GoZP8/dRNcb61Pn/JJuhygiZyHQPy8H4EyUOB9IyPJVLndCEzkz8SUL8t8uTfnvvU0pWyKODdv203fUDIZ+tIA9+4+4HZ6InKFAS3xFgJHW2o25GYxITmhSJ54GNUozfvJqPp+8mikLtjBnuVP2u6qFyn4i4SLQ/1PHAZ1zMQ6RHFUgyssdV9RkVK92NKpVltQjJ3jjy2U8Nmwayet1b7lIOAh0BDUYWGiMuRPYAJx004m1tl0OxyWSIxJKFeSZLk2Zu9yZ7bdh2376jJxBu0YV6Xx1bYoX1mw/kVAVaIIaBxwEJgKpuReOSM7zeDw0rZtAgxql+XzyasZPXsPk+ZuZs2wbHTvUokPzyir7iYSgQBNUY6CptXZJbgYjkptioiPp2L4W7RpV5PUJS1mwcgevTVjKT3Oc2X61qpRwO0QR8RPon40WKJabgYgES7lShXj2vmY81bkJZYrHsm7rPnqP/IXhnyxi74GjbocnIj6BjqAGAO8aY0YCa4Hj/p3W2u9yOjCR3OTxeGheL4ELTWk+S1rNF1PW8PO8Tcxato27OtSiffPKeCOy3vInIsEUaIL62Pd9yCn6dKOuhK2Y6Eju6lCLxEYVeW3CUhbaHYz9YgmT5myk2031qXmeyn4ibgkoQVlrdQVZ8rRypQvx3P3NmLV0G29+vYx1v++j14hfuKxJJTpdVZuihf62wpeI5LJ/2m7jjBhjNN1cwpbH46FF/XKM7tWOWxKrE+n18NPcTXQdmMR3M9eTlp63F1YWCTXZjaAeN8b0BUYAP1trj5/qIN+mhlfj7AmVCkzO8ShFgiimQCR3X1mbxMaVeO2LJSxatZMx431lvxvrY1T2EwmKbLfbMMbcADwHnIezz9NyYBfOgrGlcXbZbQ5sAv7PWvt57oZ75rTdhpyLjIwMZi7dxptfLWPX3sMAKvuJ5KCz3m7DWjsBmODbKfdKnGRUFmcliRRgATDAWvtLzoct4j6Px0PL+uVoaMrwv6RVTJi6hp/mbmLW0m3cfWUtLm+m2X4iuSXQSRJT0U65ko9llv3a+Wb7LV61k9GZZb+bGlCjUnG3QxTJczQ7T+QMVChTmOf/3Zy+dzemVNEY1mzZR88R0xn52WL2HdRNviI5SQlK5Ax5PB5aNijH6D6J3NyuOt4IDz/O3ki3QUn8MGuDZvuJ5BAlKJGzFFsgkk5X1WZEj7ZcUL00B1KPM+rz3+g1YjqrNu1xOzyRsKcEJXKOKpYtzPMPNKfP3Y0oWTSG1Zv3/ln223/omNvhiYStQJc6whhTBqgPROFMM/+T1uKT/M7j8XBxg/I0rFmWT3+yfDltLT/O3sjMJdvodFUtLmtyHhGa7SdyRgJKUMaYLsBonOSUldbiE/GJLRBJ56vrkNi4EmO/WMKSNbsY+dlvTJqzka431qd6Rc32EwlUoCOoXsAbwJPW2gNn+2bGmMuAgUB1YAcw2Fr7mjEmGhgJ3AykAUOttQP8zrsVeBFIAKYBna21O842DpHcVrFsYfp3bcGMxVt58+tlrNq0lx7Dp3NFs8rc1aEWRQpGux2iSMgL9BpURWD4OSanisB4oD/O3lL/AgYYY64A+gEGqIqzOWInY8zdvvNqA28BnYGSwGrgk7ONQyRYPB4Pl1xYnjF92nFjm2pEeDz8MGsDXQcm8ePsjaRrtp9ItgJNUJOAxHN8r8rAR9baCdbadGvtPJybf1sCnYAXrLV7fEtdDAEe8J3XEfjGWjvDWnsEeBJoaYypfo7xiARFXEwU91xThxE92lC/WikOpB5j5GeL6f3qL6zZvNft8ERCVqAlvt+AocaYa4FVwElTk6y1vf/pBXzLIf25JJIxpgRwCTAOp3SX7Hf4SqCe73FtYL7f66QaYzb7+lcHGL+I6yrFF6F/1xb8svh33vp6GXbTHp4YPo32zZ2yX+E4lf1E/AWaoFoDc4BYnAVi/Z1xncIYUxT42veaC3zNqX6HpAJxvseFsvRl7RcJGx6Ph1YXVqBRrbJ8PMny9S/r+H7mBn79bSudr3JWUNdsPxFHoGvxtc2pNzTG1AC+whkx3YmT9PD7Dk7yOeh7fChLX9Z+kbATFxNFl2vrcmnjSoydsIRla3cz4n+L+dE3269ahWJuhyjiujO5D6oszp5PdXCuXa0A3rDWrjuD12iFk5zGAk9ZazOAI8aYFJxJEr/7Dq3JXyW/ZF9f5mvEAZU4uSQoEpbOSyjCi91aMm3R77z99TLsxj30eOWvsl8hlf0kHwv0PqgmwE/AZmAmzo26VwMPG2PaWGvnZ3e+7zWqAt8CT1trX83SPQ541hizBKek1xMY7uv7CJjh2/JjFjAAWGStXRVI7CKhzuPx0OaiCjSp/VfZ77uZG/h1iVP2a9dIZT/JnwIdQb0MfAx08416ADDGjAQGA4GUALsDhXGmlg/wax8FPON7j+U4o7PXcUZZWGuXGmPu9T0vj3Pd6pYA4xYJG/5lvzFfLGH5ut0M/3QxP852yn5VVfaTfCbbHXUzGWMOAxdYa22WdgMssNYWyqX4zpl21JVwlJGRwbSFW3j7m+XsOXCUCA90aFGFju1rquwneUp2O+oGeh/UNpz7mLI6Hzjrm3dF5NQ8Hg9tGlZkbN9ErmtVFTweJv66nq6Dkvh57ibd5Cv5QqAlvnHA68aYx4DZvrbmwDBfn4jkgriYKO67ri6XNnHW9nPKfov+XNvv/PJF3Q5RJNcEOoJ6AWc1if8BW3Bm230MfAY8nTuhiUimyglFGPBgS5644yKKFS7Aig1/8Piwqbw2YQkHDx93OzyRXBHofVDHgPuNMT1xpnwfBtZYaw/nZnAi8hePx0PbhhVpUjuej35cybcz1vHtjPXMWLyVe66pTduGFfF4NNtP8o7TJihjzJXAT9ba477HWVV05khoPyiRYCoYG8X919fj0iaVGDN+CSs2/MGwjxf9OduvSjmV/SRvyG4E9S0Qj7MtxrfZHKf9oERcUKVcUQY9dDFTFmzmnW+SSV7/B48Nm8ZVLatw5xU1KRh7qu3bRMLHaROUtTbiVI9FJHR4PB7aNapEkzoJfPjDCr77dT3f/LKOXxb/zj1X16Ftwwoq+0nYCijxGGMmG2P+dpegMaa0MWbBqc4RkeApFBvFAzfUZ9jjbahVuQR7Dxxl2McL6TtqBuu37nM7PJGzkt01qDY4W12As5r5A8aYrPc81cLZZFBEQsD55YsysLuv7Pft8j/Lfle3rMIdKvtJmMnuGtRunDXxPL6v7jjbsWfKwFlRvEeuRSciZywiwkNi40o0rftX2e/rX9YxffHv3HtNHdpcpLKfhIfsrkEtxVkpAmPMFOBGa+2eYAUmIucms+x3WZPzGPuFM9tv6EcL+XH2RrrdWJ/zEoq4HaJItrIr8cVZazM3Crwqs+1Ux/odJyIhJrPsN3n+Zt6duJzl63bzyNCpXHPx+dxxhSEuRmU/CU3ZTZI4YIwp43t8EGfNvaxfme0iEsIiIjzOckl9ErmyRWXIyOCr6WvpOjCJqQu3EMii0SLBlt01qHbAH77HObajroi4p1BcNN1uasBlTZ2yn924h5c/XMCPszfQ9cb6nBevsp+EjoC228jKGBMN1AdWWWv353hUOUjbbYicWnp6BknzNvHuxGT2HzpGRISHay85n39drrKfBM85b7dhjKlmjJlmjGnmuw411/e10RjTLKcDFpHcFxHhcUZSfRPp0KIyGRkZfDltLd0GJTFNZT8JAYGuEPEqzrWmDcBdQAWcRWPHAENzJTIRCYrCcdE8eFMDhj7aGlOpOH/sP8qQDxfw9JiZbEwJ6QKJ5HGBJqhLgMettSnA9cBEa+1q4A3ggtwKTkSCp1rFYrz08CU8fOsFFI6LZunaXTz68lTe/mY5qUe0pYcEX6AJ6ggQZYwpiLOqxPe+9nhA66iI5BERER4ub3oerz2ZSIfmlUnPyGDC1DV0GzSZ6YtU9pPgCjRB/YgzWhoPpALfGGMSfW1f51JsIuKSwnHRPHhzA15+tBXVKxbjj/1HGPzBAv4zdiabt+vOEgmOQBPUA8B8nJHUVdbaQ0BjYCrwWO6EJiJuq16xOEMeacVDtzSgcFw0S9bs4uEhU3jnm+UcPnrC7fAkjzvjaebGmCJAhLV2b+6ElLM0zVwkZ+w/dIz3v0tm0pyNZGRAyaIxdLm2Lhc3KKe1/eSsnfM0cwBjTDdjzGZgD7DbGLPNGNM3RyMVkZBVpGA0D91yAUMeaUW1isXYve8IL42bzzOvzVLZT3JFoPdB9QQG4kw3vwRoBQwDehtjHs298EQk1NSo5JT9ut/cgMJxUSxevZNHXp7Cu9+q7Cc5K7uljvx1B7paaz/2a/vVGLMR6A8Mz/HIRCRkeSM8tG9emeb1Ehj3/QomzdnI+ClrmLZwC12uq0vL+ir7ybkLtMRXGph3ivYFODftikg+VLRQgZPKfrv2HWHQ+yr7Sc4INEEtA245RfttwMqcC0dEwlFm2e/BmxtQKPavst97E5M5orKfnKVAS3zPABONMc2BWb625kB74MbcCExEwos3wkOH5pVp4Sv7/Th7I59PXs3UhVu479q6tKifoLKfnJGARlDW2klAInAUZy2+m4H9QGNr7be5F56IhJu/yn6XUK1CUXbtPczA9+fxzOuz2LJDZT8J3FlttxFOdB+UiHvS0jOYNHsD73+3goOHjxPp9XBDm2rcmliDmAKBFnAkL8vuPqhst3wHXsEZLR0FJgB9Q33/JxEJHd4IDx1aVKFF/XK8NzGZn+Zu4rOk1UxZsIX7r6tL83oq+8npZVfi6wdcA7yEs6XGVThr74mInJGihQrwyG0XMvjhSzi/vFP2G/DePJ57YzZbdx50OzwJUdklqJuBO6y1A621g3Fm8V1njNFWmyJyVmpWLsHQx1rT9cb6FIyNYqHdQffBUxj3/QqOHNNsPzlZdgmqAidPIZ/nO75srkYkInmaN8LDVS2rMLZPIpc2rsSJtHT+9/Mqur80mVlLt2lLD/lTdgnKC6RlPrHWZuBci4rO7aBEJO8rVrgAj97uK/uVK8qOPYd58d25PPemyn7iCHixWBGR3FCzcgmGPt6arjfUo2BMJAtXOmW/D1T2y/f+aZ5nZ2OM/58ykUBHY8wu/4OstaNzPDIRyTe8ER6uuvh8WjYoz7sTl5M0bzOf/ryKKQs2c9919WhWN16z/fKh7BLUJqBblrYU4J4sbRmAEpSInLNihQvw2O0XcXnT8xj7xRLWb93Pi+/OpVGtstx/fV3KlSrkdogSRKdNUNbaykGMQ0TkT7WrlGTYY635buYGPvxhBfNXbGfxqp3c1K4atyTWoECU1+0QJQh0DUpEQpLXG8E1l5zPmL6JtGtUkRNp6Xz60yoefGkyc5Ztczs8CQIlKBEJacULx/D4vy5iYPeLqZxQhB1/pNL/nbn0e3M223Ydcjs8yUVKUCISFuqcX5JXHm/N/dfVJS4mkvkrttN98GQ++nElR4+n/fMLSNhRghKRsOH1RnBtq6qM7ZNI24YVOH4inY8nWbq/NJm5ySluhyc5zJXlhI0xTYBvrbVlfM+jgZE4yyulAUOttQP8jr8VeBFIAKYBna21O4IeuIiEhOJFYnjijoZc0awyY79YwoZt+/m/t+bQpHY8919fl/iSBd0OUXJAUEdQxhiPMeY+YBInr0jRDzBAVaAx0MkYc7fvnNrAW0BnoCSwGvgkiGGLSIjKLPvd5yv7zU1O4cGXJvOxyn55QrBLfP1w7q3qn6W9E/CCtXaPbz+QIcADvr6OwDfW2hnW2iPAk0BLY0z1IMUsIiHM643gulZVGdMnkTa+st9HkywPDZ7MPJX9wlqwE9RYa21DYH5mgzGmGE7pLtnvuJVAPd/j2v591tpUYLNfv4gIJYrE0OOOhgx4sCXnxRcmZXcqz781h/5vzyFlt2b7haOgJihr7bfvWckAAA95SURBVNZTNGfeGp7q15YKxPn1p3Iy/34RkT/VrVqKV55oQ5dr6xJbIJI5y1Po/tJkPvnJckxlv7ASCrP4Mv+0ifVriwMO+vXHcjL/fhGRk0R6I7i+dVXG9k2k9YUVOHYinQ9/WMlDg6cwf8V2t8OTALmeoKy1e3DW+DN+zTX5q6yX7N/n24q+EieXBEVE/qZEkRh6dmzIi91aUrFsYbbtPkS/N2fzwjtz2P5H1sKMhBpXppmfwjjgWWPMEpySXk9guK/vI2CGMaYNMAsYACyy1q5yI1ARCT/1qpViRI82fDtjHR/9uJLZy1JYuHIHt15agxvaVCNaa/uFJNdHUD7PAMuA5Tg7944HxgJYa5cC9/qe7wLq4Gw/LyISMKfsV40xfRJpdWF5jp1I54MfVvLQkCksWKmyXyjy5PXtlY0xlYH1SUlJVKhQwe1wRCRELF2zizFfLGHz9gMANK+XwH3X1qVMCc2/CqYtW7aQmJgIUMV3m9GfQmUEJSISVJllv3uvqUNsAS+zlm6j20uT+fRny/ETmu0XCpSgRCTfivRGcEMbX9nvgvIcO57GB987s/0WrtRqam5TghKRfK9k0Vh63dWI/l1bULFsIbbuOsSzb8zixXfnsmOPZvu5RQlKRMSnQfXSDH+iLfdcXZuYaF/Zb9Bk/vfzKpX9XKAEJSLiJyoyghvbVmdMn0QublCOY8fTGPf9Cqfsp00UgkoJSkTkFEoVi6XP3Y35vweaU760r+z3+iwGvKeyX7AoQYmIZOOCGmV4tWdbOl/llP1mLtnGgy9N5rOkVRw/ke52eHmaEpSIyD+IiozgpnbVGd07kZYNynH0WBrvf7eCh4dMYfEqlf1yixKUiEiAShePpa9f2e/3nQf572uzGPj+PHbuOex2eHmOEpSIyBnKLPt1uqo2BaK9/PrbVrq9lMTnk1er7JeDlKBERM5CVGQEN7erzpjeibSs75T93puYzMNDpvDbqp1uh5cnKEGJiJyD0sVj6dupMf3+3ZzypQvy+86D/Oe1mQx6fx679qrsdy6UoEREcsBFxin73X1lLaKjvMz4bSvdBiUxXmW/s6YEJSKSQ6IivdySWIMxfdrRvF4CR46l8e7EZB4dOoXfVqvsd6aUoEREcliZ4nE81bkJ/e5vTrlSBdm8/SD/GTuTwePms3ufyn6BUoISEcklF9Usw8hebbmrg1P2m774d7oNSmLC1DWcSFPZ758oQYmI5KKoSC+3XlqDMb2dst/ho2m8/c1yHnl5KkvWqOyXHSUoEZEgKFPCKfs9d38zEkoVZPP2Azw9ZiaDP1DZ73SUoEREgqhhzbKM7NmWjh1qOmW/RSr7nY4SlIhIkEVHebntUsPo3u1oVjf+pLLf0jW73A4vZChBiYi4pGyJOJ6+pynP3teMhJJO2e+pMb8y5IMFKvsBkW4HICKS3zWqVZb61UrxxdQ1fPbzKqYt2sLc5BTuuKImV19chUhv/hxL5M+fWkQkxERHebn9MsOo3u1oWieew0dP8NbXy3hs6FSWrc2fZT8lKBGREBJfsiD/ubcpz3RpSnzJODamHODJ0b/y8kcL2LP/iNvhBZVKfCIiIahx7XgaVC/N+Clr+DxpFVMXbGHucl/Zr2UVvPmg7Jf3f0IRkTAVHeXlX5c7Zb8mteNJPXKCN79axmPDprF83W63w8t1SlAiIiEuvmRB/tulKf/t0pSyJeLYsG0/fUfNyPNlPyUoEZEw0aR2PKN6t+OOyw1RkRFMXbCFroOS+Hr6WtLy4E2+SlAiImGkQJSXf11Rk1G92tGoVllSj5zgjTxa9lOCEhEJQwmlCvLsfc34771NKeNX9hv28UL2HMgbZT8lKBGRMNakTjyje7fj9sucst/k+ZvpOjCJb35ZF/ZlPyUoEZEwVyDKy53tTy77vf7lUh5/ZRrJ68O37KcEJSKSRySUKsgzXZryn3uaUKZ4LOu37qfPyBm88slC9h446nZ4Z0wJSkQkD/F4PDStm8Co3u247bIaRHojSJq3ma4Df+bbGeFV9lOCEhHJg2KiI+nYvhajerelYc0yHDpygtcmLOWJV6azYv0fbocXECUoEZE8rFypQjx7XzOe6tyE0sVjWbd1H71H/hIWZT8lKBGRPM7j8dC8XgKje7fj1kv9yn6Dkpg4Yx1p6Rluh3hKSlAiIvlETHQkd3WoxahebbnIlOHQ4eOMnbCUJ16ZxsoNoVf2U4ISEclnypUuxHP3N+Opzo2dst/v++j16i+M+HQR+w6GTtlPCUpEJB9yyn7lGN2rHbckVifS6+GnuZt4YGAS381cHxJlPyUoEZF8LKZAJHdfWZuRvdpxQY3SHDp8nDHjl9Bj+DRWbnS37KcEJSIilC9diOf/3Zy+nRpTqmgMa7fso9cId8t+SlAiIgI4Zb+W9csxpk8iN7f7q+zXdWAS37tQ9guLLd+NMQ2AsUB9YB1wr7V2nrtRiYjkTTEFIul0VW3aNarI6xOWsnj1TkaPX8KkORvpdlMDalQqHpQ4Qn4EZYyJBr4CPgWKAS8Ak4wxRVwNTEQkj6tYtjDPP9Ccvnc7Zb81W/bRc8R0Xv3f4qCU/UI+QQFtgChr7SvW2uPW2k+A5cBt7oYlIpL3eTweWjYox+g+idzUthreCI8zkhqUxA+zNuRq2S8cElRtYEWWtpVAPRdiERHJl2ILRNL56jqM6NGWC6qX5kDqcUZ9/htPjZ7BseNpufKe4ZCgCgGpWdpSgTgXYhERydcyy3597m5EyaIxrP19H4ePnsiV9wqHSRKHgNgsbXHAQRdiERHJ9zweDxc3KE/TOvEcPppGkYLRufI+4TCCSgZMlraavnYREXFJVKQ315IThMcIagrgMcY8DowEbsKZbj7B1ahERCRXhfwIylp7DOiAk5j+AJ4GrrfW7nQ1MBERyVXhMILCWrsMuNjtOEREJHhCfgQlIiL5kxKUiIiEJCUoEREJSWFxDeoceQFSUlLcjkNERLLw+93szdqXHxJUAsCdd97pdhwiInJ6CcBa/4b8kKDmAZcA24DcWTBKRETOlhcnOf1tCyVPRob7+86LiIhkpUkSIiISkpSgREQkJClBiYhISFKCEhGRkKQEJSIiIUkJSkREQpISlIiIhCQlKBERCUn5YSWJs2aMaQCMxdnBdx1wr7X2b3c7y6kZY+4FXgOO+jV3t9a+51JIIc8Y0wT41lpbxvc8Gmcn6ZtxVkIZaq0d4GKIIe0U/34FgAPAMb/DZlprL3cjvlBkjLkMGAhUB3YAg621r4XCZ08J6jR8/3G+Al4BWuHs6DvJGHOetXa/q8GFj4uAl621fd0OJNQZYzxAF2BIlq5+gAGqAkWBH4wxv1tr3w9yiCEtm3+/esAf1tr44EcV+owxFYHxQCec33cNgR+NMRuANrj82VOJ7/TaAFHW2lestcettZ8Ay4Hb3A0rrDQEFrsdRJjoB3QD+mdp7wS8YK3dY63dgPML+IEgxxYOTvfvp89g9ioDH1lrJ1hr030VoqlAS0Lgs6cEdXq1gRVZ2lbi/EUm/8AY48Upjd5ljNlqjFljjOnr+0tX/m6stbYhMD+zwRhTDGcRzWS/4/QZPLW//fv5XASUMcYsMcZsN8Z8Zowp70J8Icla+4u1tmvmc2NMCZzFtRcRAp89JajTKwSkZmlLBeJciCUclcb5ZfEeUAWnjt3N9yVZWGu3nqK5kO+7/+dQn8FTOM2/H8Ah4FcgEadcdRiYEKy4wokxpijwNTAHWOBrdvWzp2tQp3cIiM3SFgccdCGWsGOtTQFa+zUtNsa8inMtb7Q7UYWdQ77v/p9DfQbPgLX2Cf/nxpgngJ3GmIrW2s0uhRVyjDE1cK5BJQN38tdnztXPnkZQp5eM8xeXv5qcPOSV0zDG1DHG9MvSHA0ccSOecGSt3QOkcPLnUJ/BM2CMed4YU8uvKdr3XZ9DH2NMK5xR05fAzdbaI6Hy2dMI6vSmAB5jzOM4Uy1vwrmmovJAYPYCPYwxW4C3gAuBR4CHXI0q/IwDnjXGLMEp+fUEhrsbUlipDzQyxtzhez4cmGit3eliTCHDGFMV+BZ42lr7apZu1z97GkGdhrX2GNABJzH9ATwNXK8PdmCstb8D1+LM+tmPM5X1/6y1n7saWPh5BliGM4N0Hs6/41hXIwovXYA9wBpgA879UHe5GVCI6Q4UBgYYYw76fQ0iBD572lFXRERCkkZQIiISkpSgREQkJClBiYhISFKCEhGRkKQEJSIiIUkJSkREQpJu1BXJIcaYd3FWgD6dfjgrRU8BCltrg7JsjG/h3l+Bu621q7I5LgKYDdxlrbXBiE0kOxpBieScR3FWgE7A2a4FoIlf2xBgpu/xoVOcn1seAX7LLjkBWGvTgefRjcASInSjrkguMMbUBZYCVXx76bgVRwywCWhnrV0W4DlrgS7W2qm5GZvIP1GJTySIjDFt8CvxGWMygH8BT+IszDkf6Aj0wlmSZz/wpLV2nO/8wsDLONuXZACTgUez2W7idmCvf3IyxvwX+DfOligrgKestd/7nTMBZzQ4NQd+ZJGzphKfiPsGAo8BzYBKwEKcxNQY+AJ4zRiTuTfU6ziJ7Aqc7UwycLboPt0fm1cBP2Q+Mcbc4HuvjjirU08EPjPGFPE75wfg0mxeUyQolKBE3DfKWjvFWrsYZ2XpgzijGgsMxdmTp4ox5nycEdEd1tp5vlHRXTjbdrc/zWs3wlnsM1Nl4Ciw0Vd6fB64ETjud0wyzurVNXPkpxM5S/oLScR9a/wepwIbrLWZF4cz9y0qAJzne2yNOWmrsjicUdW3p3jtssAuv+cf4Mw0XGeMWYCzg+o71trDfsfs9n0vc4Y/h0iO0ghKxH3HszxPP81xkb5jLwQu8PuqAbxzmnPSAU/mE992MQ1xRlwzgc7AEt+kjkyZvxfSAv4JRHKBEpRI+FgBRAEFrbVrrLVrgG3AYJwkdSopOJMhADDG3Ag8YK2dZK19FGfkdQC40u+c0n7nirhGJT6RMGGttcaYr4H3jTHdgZ3ACziTK1ae5rQFQAO/515gsDFmO86MwWZAvO9xpgb8tcmfiGs0ghIJL51wksmXOLucFgUus9buPc3xE3Fm+wFgrf0MeBZn1LUK6A88ZK2d7HdOK+AHa61KfOIq3agrkocZY+Jwtjpvb61dGMDxEcBGnJmCv+RyeCLZ0ghKJA+z1qbijJa6B3jKdcA6JScJBUpQInnfMKC+yTI3PSvf6OlpoGtQohL5ByrxiYhISNIISkREQpISlIiIhCQlKBERCUlKUCIiEpKUoEREJCT9Pz/ZkiJo7GgeAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plot_position(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3deZicVbXv8W9V9ZB0OglJSEwCBAnCAmQMyMEoEREHOBEMcBEMHDnguQ4oUePhBEWZDCKg4oACR0REpsskEBTDEB4ZBKJIkAQWhsEgZB5I0mN1Vd0/3upOpenuVCf1DtX5fZ6nn+7a71tVK6Holb332nunCoUCIiIiSZOOOwAREZGeKEGJiEgiKUGJiEgiKUGJiEgi1cQdQNjMrB54H7AUyMUcjoiIbC4DjAPmu3tb6YUBn6AIktNjcQchIiJ9Ohx4vLRhe0hQSwFuuukmxo4dG3csIiJSYtmyZUyfPh2Kv6tLbQ8JKgcwduxYdt5557hjERGRnr1jCkZFEiIikkhKUCIikkhKUCIikkhKUCIikkiJK5Iws68BH3L3T5W0GfALYBKwAbja3WfHFKKIiEQgMT0oM2s0s8uBH/Rw+SbgIWAkcCRwtpkdG2V8IiISrST1oO4HVgLXEKwqLmXF7ymgUPxqjS40iUshn4N8vvg9R6HkZwp5CoU8FAoU8nko/kyh0NVOIfi4FEp+3vSdTdeCd+tqCx6WHEVT/LnQda3HaHto2srjbHQMjlSJmuGjqRu9SzivHcqr9sDM6gh6QN0V3H05cIq7v2VmF/DOBHUxMBu4iGBbjB+5+9ww45XyFQp58q3N5FvWk2veQK5lA/nWJvJtLRTam8m3tZBvayafbaOQbaWQbSdf/F7oyFLIZYPvHe0U8h0Uch2QywXfe84EIpIgE75yLTXDRlX8daPsQU0G5vXQngNq3P2tPp5bAGYC1xL0pu4xs4Xufl3lw5Tu8m3NZNcsJbvmLbLrVpBbv5qODWvo2LCG3IbV5JrXB72XUKQgkyGVzkA6QyqdhlQ6eJxKk0qloKst+E4qFbST6noMqeK34qh213WK11Ilb1m8Bl3XN7V3Pep6/I54y2jqpXEr7hGJV+0OY8g07hDKa0eWoNz9Ubbi/zgzOwT4mruPLzYtMLPLgC8BSlAVVMjnyK5+k7alr9K29BXal79Gds1b5Jre3uJzU/UNZAYPJdMwjPTgRjKDGknVDyZd3xB81Q0mXTeIVG09qZq64OeaOlI1tcH3TM2m75maIBllaoJEJCLbpSTNQfVmF6DOzFLu3jne0wFkY4xpQCjkc7Qte42WV5+j5bUFtC19lUL2nVN7qZo6akaMpXbkOGpHjKVm2CgyQ0dSM3QUNUNHkhkynFSmNoY/gYgMZNWQoJ4gmHe60MwuAnYDvgFcFWtUVSqfbaPZn6HpH/Npee158i0bNrteM3w0dWMnUj/uPdSP3Y26HXcmM2wUqVRiCj5FZDuR+ATl7ivM7GjgcuBsYC3BXNRPYw2sihQKBdqXvsL6BQ/TtPBx8m3NXddqho9h8MQDaZh4IIMm7E2mYViMkYqIbJK4BOXuF/TQ9hTBWSHSD4V8jg3Pz2P9/N/TvuKfXe314/egcd/DGTzxIGpHjtu8QEBEJCESl6Bk2xUKBZr8KdY+ejPZ1UFxZLphGEP3ncLQAz5C3ZgJMUcoIrJlSlADTMtrz7Nm3k20LV0MQM2IsYyY8mka936/ChlEpKooQQ0Q+bYWVj1wLRtf+BMAmSE7MOLwkxh64EeCsm0RkSqj31wDQNvSV1l+9w/oWLuMVG09O3zgRIa/7xjSdYPiDk1EZKspQVWxQqHA+r/8gdUP3wC5DurG7MqY42dSN2qnuEMTEdlmSlBVKt/ewop7fkLzy88AMOzgTzDyqM+SrqmLOTIRkcpQgqpC+Wwby277Hq1LFpKub2DHqV+ica/3xx2WiEhFKUFVmXxHO8tvv5TWJQvJNI5k/GkXUjty/JafKCJSZbR/TRUpdGRZfsfltLz2PJkhOzDu1AuUnERkwFKCqhKFXAfL7/4BLa88S7phGOOmn69iCBEZ0JSgqkChkGfFvT+h+eX5pAc1Mu6U71A3WrtBiMjApgRVBdb/5QGaFj1Bqr6Bcad8m/qxu8UdkohI6JSgEq591b9Y88iNAIyZehb1498Tc0QiItFQgkqwQq6DFff8hEJHO437H8GQvQ6LOyQRkcgoQSXY2sfvoH3ZK9QMH82OHz0j7nBERCKlBJVQrW++zLon7gRSjP7kV0gPGhJ3SCIikVKCSqB8eysr7vkxFPIMP+xYBu/63rhDEhGJnBJUAq155EY61i6jbswERn7olLjDERGJhRJUwmTXLmP9s3MhlWb0sTNI1eiQQRHZPilBJcy6J+6CQp7G/T5E/bveHXc4IiKxUYJKkOy6FWz4+6OQSjPiAyfEHY6ISKyUoBJk3ZN3QT5H476HUztyXNzhiIjESgkqITreXsmGBfMglWYH9Z5ERJSgkmLdk3dDvoPGfT6gXcpFREjQgYVmNgOYAYwCHJjp7o8Vr00ArgMOA1YAX3H338cVa6V1rF/N+gUPAyl2+OCJcYcjIpIIiehBmdnxwDnAVGAE8AtgjpmNLt5yK/A8QfL6L+BWM5sYR6xhWPfnuyHXwZB9JlO3485xhyMikgiJSFDAOOASd1/k7nl3vx7IAfuZ2Z7AIcB33L3d3R8B7gXOjDHeiunYsIYNf3sIgBEfUO9JRKRTZEN8ZlYHjOzhUsHdr+p27xSgEVgIvB9Y4u5NJbe8BBwaVqxRenv+/RRyWYbsdRh1Y3QIoYhIpyh7UJOBpT18vVl6k5ntC9wGnOfuywkSVXO312oGGsIOOGyFQp6NL/wJgOH/dmzM0YiIJEtkPSh3fxRI9XWPmU0FbgQudffLis1NwOButzYAGysdY9Ral7xIbsMaanYYQ/1Oe8YdjohIoiRlDqqziu8W4HPu/v2SS4uACWZWmqT2KrZXtY0LHwegcZ8Pkkr1mbtFRLY7iUhQZnYScAlwlLvfWXrN3R1YAMw2s3oz+zBwHHBz9JFWTiHXQdNLfwag8b0fjDkaEZHkSco6qFlAPfCwmZW2n+zuc4ATgGsJ1kCtAs509xcij7KCWl5bQL5lA7Wjd6FuzK5xhyMikjiJSFDuPmkL198Ajo4onEiUDu+JiMg7JWKIb3uTz7bR9PIzgIb3RER6owQVg+bFf6XQ3kr9+D2oHTE27nBERBJJCSoGXcN76j2JiPRKCSpi+dYmWhY/C6QYsvfkuMMREUksJaiINfnTFHJZBu36XmqG9rTzk4iIgBJU5DYu0vCeiEg5lKAilGt6m5bX/g7pGobsdVjc4YiIJJoSVISa/jEfCnkaJh5AZvDQuMMREUk0JagItS55EYDBEw+MORIRkeRTgopQ65Jgf9tBE/aJORIRkeRTgopIx9sr6Xh7BelBQ6gbvUvc4YiIJJ4SVERa3giG9wbtvBepdCbmaEREkk8JKiIa3hMR6R8lqIgoQYmI9I8SVARyTW+TXf0mqdp66sdOjDscEZGqoAQVgZY3ir2nnfYklUnEEVwiIomnBBUBDe+JiPSfElQEOhfoKkGJiJRPCSpkudYm2pe/Duka6sfvEXc4IiJVQwkqZG3/egkoUD/+PaRr6+MOR0SkaihBhaylOP80WMN7IiL9UlZJmZntBxwNHAKMAXLAMmA+MMfdF4cWYZXT/JOIyNbpswdlZlPMbB7wF+CTwFrgyeLjVuBUYJGZPWhmU8IOttrks220LV0MqTSDdra4wxERqSq99qDM7FfAe4GfAdPcfV0v9w0DTgGuNLPn3f30MAKtRm1vvgz5HHVjJ5Kub4g7HBGRqtLXEN/97n7Gll7A3dcD1wDXmNmJFYtsAGjR+icRka3Wa4Jy9zv7+2LufsfWBmJmM4AZwCjAgZnu/ljx2sHAj4D9gfXAL4GL3b2wte8Xhc4FuoN3UYISEemvcoskGoD/AX7r7v8ws2uA6cAzwHR3X7otQZjZ8cA5wEeBl4DPAnPM7D1AE3A/8F3gw8BE4I8ERRrXbsv7hqmQywZDfMCgXfaKORoRkepTbpn5jwkKIurM7DiCBPI1oAX4SQXiGAdc4u6L3D3v7tcTVAruB+wC/Nndf+buOXf/B/A74IMVeN/QtK98g0JHO7Ujx5MZMjzucEREqk65O5ceB0x194Vmdi7woLv/r5k9Afy5nBcwszpgZA+XCu5+Vbd7pwCNwEJ3Xw5M6/Y6R5Pg3hNA+6p/AVA3ZkLMkYiIVKdyE9RgYLmZpYGPAxcU2wsEPZ1yTAbm9dCeK43DzPYFbgPOKyYnSq7VAzcDzcDVZb5vLLIr3wCgdsedY45ERKQ6lZug5hPMQa0ERgB3m9l44GLgqXJewN0fBVJ93WNmU4EbgUvd/bJu18YCdwJ54Ch3bykz9li0rwoSVN2Ou8QciYhIdSp3DurLBD2gGcCX3P0t4FzAgLMrEUixiu8W4HPu/v1u1/YhSJKLCZLT2kq8Z5iynUN8o5WgRES2Rl8LdScDTxWLFhYBB3a75ZvuvqESQZjZScAlwJHu/nS3ayOAucCt7v6NSrxf2PId7WTXLodUmtqR4+MOR0SkKvU1xHcFYGb2GEGCmFu6516lklPRLKAeeNhssy2BTiYoK98J+KKZfaHk2n3ufkoFY6iY7Oq3oJCndtR4UjW1cYcjIlKV+lqoO7m4jdFHCNYnfd3MMsCDBAnr4UoNtbn7pC3cUolS9shkV3UWSGh4T0Rka/VZJFHcxuju4hdmthvwMYKezdVm9gpByfl5YQdaTdpXFuefVMEnIrLVyq3iA8DdX2PTvntp4H0EvSsp0VXBpwIJEZGtVnaCMrMjCHY3734sbHMlAxoINMQnIrLtyt2L70rgK8ASgnOgShWAH1Y4rqpV6MiSXbMsqOAbpQo+EZGtVW4P6j+AM9z9hjCDGQiya4oVfCPHka6pizscEZGqVe5C3WaCnctlCzr34NMWRyIi26bcBPVd4IpiFZ/0oX2ltjgSEamEcof4XiTY6WFxt4W0ALh7ppJBVbOuAglV8ImIbJNyE9S1BJvCXo+q9vrUdcyGhvhERLZJuQlqF+Bod381zGCqXSGXJbtmKZCidtROcYcjIlLVyp2DehCYEmYgA0F2zVLI56gZ8S7Std2Xi4mISH+U24N6GviZmZ1AcORFtvSiu59T6cCqkYb3REQqp9wE9VGC85gaeeexG4WKRlTFsitVYi4iUillJSh3/3DYgQwE7auWANqDT0SkEnqdgzKz881scLkvZGZDzeyiyoRVnTYN8SlBiYhsq756UG8DC83sDuAud3+q+w1mlgIOAU4Fjgd+FEqUVaCQ6yC7eimAKvhERCqgrwMLrywmp3OAuWbWQbBgdxWQAkYT7G6eAn4NfMDdl4QecUJl1y6DfAc1O4whXTco7nBERKrelg4s/BdwtpmdCxwBHAy8C8gTVPZdCMxz97aQ40y8rjOgNLwnIlIR5RZJNAH3F7+kB6rgExGprHIX6soW6BRdEZHKUoKqEJ2iKyJSWUpQFVAoFMiuXQFA7chxMUcjIjIwlJWgzOzdIcdR1fJtzRSyraTqBpGub4g7HBGRAaHcrY4Wm9lTwG+B2919dYgxVZ3c+uCvo2boKFKpVMzRiIgMDOUmqInAZ4AvAj82s7nATcA97t5SiUDMbAYwAxgFODDT3R/rdk8twblU97n7BZV430ro2FBMUMN2jDkSEZGBo6whPndf4u6XuvsBBGuhFgDnAsvN7AYzO2pbgjCz4wkWBE8FRgC/AOaY2ehut36Xd25WG7uO9asAyAwdFXMkIiIDx9YUSfwLeAV4laAHdiDwGzNzM3v/VsYxDrjE3Re5e97drwdywH6dN5jZEQS7qv9xK98jNB2dQ3zDRsYciYjIwFHWEJ+ZDQGOA04GPgYsB24GznP3hWaWBn4O3AZM6OU16oCefoMX3P2qbvdOITjaY2Hx8Qjgf4FpBL2oRMlt2DQHJSIilVHuHNQKoB24C/iEuz9aetHd88V5qcP7eI3JwLwe2nOlcZjZvgSJ7jx3X15svhr4ubu/YGZlhhydTT0ozUGJiFRKuQnqdODenvbcM7Mx7r7C3e8iSGA9Kia1PkvczGwqcCNwqbtfVmw7HdgRuLLMWCO3qUhCPSgRkUopdw7qVmBY90Yzm0AwF7XNilV8twCfc/fvl1w6BTgUWGtm64B/B2aZ2ZxKvG8ldPagVCQhIlI5vfagzOwUgjkfCHo+vzSz7j2oXYE12xqEmZ0EXAIc6e5Pl15z9493u/d3wHNJKTPPtzVTaG8hVVtPetCQuMMRERkw+hrie5Cgaq5zWK6l+NWpQHDkxq8rEMcsoB54uNsc08nunpieUk86S8xrhmmRrohIJfV1YOEq4AwAM3sduNzdm8MIwt0n9ePeT4URw9bqWK8KPhGRMPQ1xHcM8KC7Z4H5wBG9VdC5++/DCS/5OgskMiqQEBGpqL6G+OYAYwlKzPsaZisAmUoGVU1y64MpOPWgREQqq68hvnRPP8vmNs1BaQ2UiEgllZ14zOwMMzux5PH/M7PTwgmrenRoFwkRkVCUex7Ut4Ar2Hwo7+/AlWb21TACqxaagxIRCUe5PajPE5R839bZ4O4XA6cSHJGx3cqpik9EJBTlJqgRwD97aH8FeFflwqku+bYW8m3NpGrqSA9ujDscEZEBpdwE9RTB9kKlm7pmgJkEJejbpdI9+LRIV0SkssrdLPYbwMPAEjN7nqC0fL/i848OKbbE0x58IiLhKfdE3QWAAbOBxcCLBOcy7eHufwsvvGQr3eZIREQqq9weFO6+2sz+CLxBkNhecvcNoUVWBXRQoYhIeMo9UbcRuA44EcgSbCBbY2YPAie4e1N4ISbXpoMKlaBERCqt3CKJHxLMOb0fGAwMKv48Hvh+H88b0LrWQKkHJSJSceUO8R0PTHP3Z0ranjGzs4A7gC9XPLIqoB6UiEh4yu1BpYFVPbSvAbbbBUBdc1Dah09EpOLKTVB/Ai4ws7rOBjOrB84HHgsjsKTLt7eQb20ilaklPXho3OGIiAw4/VkH9Tjwhpk9V2w7AGgFPhFGYEnXsSE4ZiOjRboiIqEodx3UYmBv4GKC7Y1eAM4D9nb3l8ILL7m0BkpEJFz9WQe1FvhZiLFUFW0SKyISrr6OfJ9PsKXRFrn7oRWLqEp0DvGpByUiEo4tHfkuvegc4tMaKBGRcPR15PuFUQZSbTatgVKJuYhIGMqegzKzk4D/BvYAJgFfApa5+xUhxZZo2odPRCRc5R75fjrwc+AuoHMt1EvAd8xsVjihJVvpWVAiIlJ55S7UnQl80d2/B+QA3P2XwH8SHAe/Xcln28i3bIRMDekGLdIVEQlDuUN8uwN/6aH9OWBsJQIxsxnADGAU4MBMd3+seG0o8FPgWILKwjuAL7t7thLv3V8dJSXmqVS5OV5ERPqj3N+uDhzVQ/tJBEN928TMjgfOAaYCI4BfAHPMbHTxll8BOwDvJlgwfAjBfFgschreExEJXbk9qG8Cd5jZIcXnfMHM3kOQUE6sQBzjgEvcfVHx8fVm9gNgPzN7ETgO2Mnd1wPrzew4IFOB990qHVqkKyISurISlLv/wcwOJei1vAB8lODY98Pc/dlyXqO40ezIHi4V3P2qbvdOIdglfSFwMLAEmG5mZwO1wG+Bb5fzvmHoWgOlHpSISGj62kniGOABd88DuPtC4PRteK/JwLwe2nOlcZjZvsBtwHnuvtzMRhIM7e0L7A+MAe4DNgCXbEM8W61DJeYiIqHrqwd1D7DazG4BbnD35/q4d4vc/VGCo+J7ZWZTgRuBS939smJzG8Fw3kx33whsNLMfAl8kpgSV27gOgMzQEXG8vYjIdqGvIomdCBLAvwHPmtnfzewbZjY+jECKVXy3AJ9z99Jj5DuLMHYoaSt7gXEY8q0bAcjoHCgRkdD0mqDcfYW7/8TdJwMTgZuB04B/mtlcMzvVzAZXIojiLhWXAEe5+53d4vg7QYn7j8xsiJntCnyNIJnFIteyAYD0oO32MGERkdCVex7U6+7+PXc/gOCgwmcIKvuWm9n1FYhjFlAPPGxmG0u+phavH0NwOOKrBMnqHuDKCrzvVsm3qAclIhK2fg+VufsiM7uCYG3UVwl6Vf+5LUG4+6QtXF8JfGZb3qNSCoVCV4JKD1YPSkQkLP3ZLHY48CmCxbkfIejN3AQcH05oyVToaKeQy5KqqSNdWx93OCIiA1afCaokKf0fgp0k1gO3Ahe4+/zww0uevOafREQi0dc6qPsJekp5gnVHJxKsi+qIKLZEyml4T0QkEn31oBqBs4Dbi1sMCZt6UCqQEBEJV18n6n4oykCqRa5VPSgRkSjorIh+yjcXe1CagxIRCZUSVD/l1YMSEYmEElQ/5bRIV0QkEkpQ/aRFuiIi0VCC6icVSYiIREMJqp+6ysxVJCEiEiolqH7atFBXc1AiImFSguon7WQuIhINJah+Upm5iEg0lKD6IZ9to9DRTipTS6qmLu5wREQGNCWofsiXzD+lUqmYoxERGdiUoPqh66h3De+JiIROCaofOuefMkpQIiKhU4Lqh5wOKxQRiYwSVD+oxFxEJDpKUP2gffhERKKjBNUPOc1BiYhERgmqHzoPK9QclIhI+JSg+qFrJ/MGzUGJiIStJu4AOpnZDGAGMApwYKa7P1a8ZsAvgEnABuBqd58ddYxdRRLqQYmIhC4RPSgzOx44B5gKjCBIRnPMbHTxlpuAh4CRwJHA2WZ2bNRx5ls7F+qqByUiErZEJChgHHCJuy9y97y7Xw/kgP2K1634PQUUil+tUQeZa1aRhIhIVCIb4jOzOoIeUHcFd7+q271TgEZgYbHpYmA2cBGQAX7k7nNDDLdHm3YyVw9KRCRsUfagJgNLe/h6s/QmM9sXuA04z92XF5sLwEyCpHUgcLyZnRlR3IB2MhcRiVpkPSh3f5RgiK5XZjYVuBG41N0vK7YdAnzN3ccXb1tgZpcBXwKuCy/izZUu0tVO5iIi4UvKHFRnFd8twOfc/fsll3YB6sysNCt0ANko49NO5iIi0UpEmbmZnQRcAhzp7k93u/wEwbzThWZ2EbAb8A3gKiK0aSdzzT+JiEQhEQkKmAXUAw8HS566nOzuc8zsaOBy4GxgLXAt8NMoA+wa4tMaKBGRSCQiQbn7pC1cfwo4PKJwepRrUYm5iEiUEjMHlXR5zUGJiERKCapMOc1BiYhESgmqTJqDEhGJlhJUmTaVmasHJSISBSWoMuVVJCEiEiklqDJpJ3MRkWgpQZVJZeYiItFSgiqTiiRERKKlBFWGzp3MydSQqq2POxwRke2CElQZNhVIDNVO5iIiEVGCKsOmgwo1vCciEhUlqDJ0FUho/klEJDJKUGUoPaxQRESioQRVhs5dJLQPn4hIdJSgyqA5KBGR6ClBlaFrH75B6kGJiERFCaoM2odPRCR6SlBl0E7mIiKbmzlzJrNnz+56nMvlmDx5Mk8//XTF3iMRR74nXb5VPSgRicfSW2fT8sqzkb3f4N0nMe7kb23xvmnTpjFr1ixmzZpFJpPhiSeeYNCgQRx66KEVi0U9qDJoHz4Rkc1NnjyZdDrd1WO67777+OQnP1nR3XbUgyqDdjIXkbiU05uJQzqd5thjj+W+++7joIMO4qGHHuKuu+6q7HtU9NUGqLzmoERE3mHatGk89NBDzJ07lz333JPddtutoq+vBLUF2slcRKRnu+++O7vuuitXXnklxx13XMVfXwlqC/KtTUCwD592MhcR2dy0adNYuXIlxxxzTMVfOzFzUGb2P8BZwEhgETDT3R8rXpsAXAccBqwAvuLuv48iLu3DJyLSu+nTpzN9+vRQXjsRPSgzOxH4MnAkMBT4FXCPmWWKt9wKPA+MAv4LuNXMJkYRm/bhExGJRyISFHAnsLe7LwYGEfSi1gB5M9sTOAT4jru3u/sjwL3AmVEEph6UiEg8IhviM7M6gsTTXcHdlwMbzewTwP1AB3CSuxfMbB9gibs3lTznJaByq8H6kGvVPnwiInGIcg5qMjCvh/ZcSRzzCHpQnwZuM7ODgEagudtzmoGGkOLcTP3Y3akZtiMNe0yK4u1ERKQosgTl7o8CfZbBuXtb8cffmtkXgKOB14DB3W5tADZWOsae1I/djQlfuSaKtxIRkRKJmIMys3PM7BfdmuuBdQQVfRPMrDRJ7VVsFxGRASopZeaPA+eb2c3An4EzgAnAve6+yswWALPN7FyCocLjgPfHFq2IiIQuET0od38S+BzwS2AVcApwlLuvKt5yArA3wRqoXwJnuvsLccQqIiLRSEoPCne/Bbill2tvEMxHiYjIdiIRPSgREZHulKBERCSRlKBERCSREjMHFaIMwLJly+KOQ0REuin53Zzpfm17SFDjgNB22xURkYoYB7xS2rA9JKj5wOHAUoJtlUREJDkyBMlpfvcLqUKhEH04IiIiW6AiCRERSSQlKBERSSQlKBERSSQlKBERSSQlKBERSSQlKBERSSQlKBERSSQlKBERSaTtYSeJrWZmBwBXA/sDrwJnuPs7VjtLz8zsDOAaoK2k+Sx3vyGmkBLPzA4F5rj7mOLjOuBnwIkEO6H80N2/F2OIidbD3189sAFoL7ntSXf/WBzxJZGZfRS4FNiD4FDYy939miR89pSgelH8j3MPcCUwheBU37lmtqu7r481uOoxCfiBu8+KO5CkM7MUcCZwRbdLFwIG7A4MBx4wszfd/TcRh5hoffz97Qescfex0UeVfGa2C3An8FmC33cHA380s9eBI4j5s6chvt4dAdS6+5XunnX3W4GFwKfjDauqHAw8F3cQVeJC4IvAd7u1fxaY7e5r3f11gl/An484tmrQ29+fPoN9ezdws7vf7e754gjRo8AHSMBnTwmqd/sAL3Zre4ngX2SyBWaWIRgaPc3M3jKzxWY2q/gvXXmnq939YOAvnQ1mtgPBJpqLSu7TZ7Bn7/j7K5oEjDGz581suZndbmY7xRBfIrn7Y+7+hc7HZjaSYHPtv5GAz54SVO8ageZubc1AQwyxVKPRBL8sbgB2IxjH/mLxS7px97d6aG4sfi/9HOoz2INe/v4AmoAngI8QDFe1AFDHO0kAAASJSURBVHdHFVc1MbPhwL3A08Bfi82xfvY0B9W7JmBwt7YGYGMMsVQdd18GfKik6Tkz+ynBXN7P44mq6jQVv5d+DvUZ7Ad3/3rpYzP7OrDSzHZx9zdiCitxzGxPgjmoRcB0Nn3mYv3sqQfVu0UE/+IqtRebd3mlF2b2XjO7sFtzHdAaRzzVyN3XAsvY/HOoz2A/mNlFZrZ3SVNd8bs+h0VmNoWg1/Q74ER3b03KZ089qN7NA1Jm9jWCUssTCOZUNDxQnnXATDP7F3AdcBBwNvDlWKOqPjcC55vZ8wRDft8AfhxvSFVlf+AQM/tM8fGPgfvdfWWMMSWGme0OzAG+5e4/7XY59s+eelC9cPd24GiCxLQG+BbwKX2wy+PubwLHElT9rCcoZb3Y3e+INbDq8x3gBYIK0vkEf49XxxpRdTkTWAssBl4nWA91WpwBJcxZwFDge2a2seTr+yTgs6cTdUVEJJHUgxIRkURSghIRkURSghIRkURSghIRkURSghIRkURSghIRkUTSQl2RCjGzXxPsAN2bCwl2ip4HDHX3SLaNKW7c+wTwH+7+ch/3pYGngNPc3aOITaQv6kGJVM4Mgh2gxxEc1wJwaEnbFcCTxZ+benh+WM4GFvSVnADcPQ9chBYCS0Jooa5ICMxsX+DvwG7Fs3TiimMQsAQ40t1fKPM5rwBnuvujYcYmsiUa4hOJkJkdQckQn5kVgFOAcwk25vwLcCrw3wRb8qwHznX3G4vPHwr8gOD4kgLwCDCjj+MmTgbWlSYnM/s28H8JjkR5Efimu/+h5Dl3E/QGH63AH1lkq2mITyR+lwJfBQ4DJgDPEiSm9wF3AdeYWefZUNcSJLKPExxnUiA4oru3f2z+O/BA5wMzm1Z8r1MJdqe+H7jdzIaVPOcB4Kg+XlMkEkpQIvG7yt3nuftzBDtLbyTo1TjwQ4IzeXYzs4kEPaLPuPv8Yq/oNIJjuz/Ry2sfQrDZZ6d3A23AP4tDjxcBxwPZknsWEexevVdF/nQiW0n/QhKJ3+KSn5uB1929c3K489yiemDX4s9uttlRZQ0Evao5Pbz2u4BVJY9/S1Bp+KqZ/ZXgBNXr3b2l5J7Vxe9j+vnnEKko9aBE4pft9jjfy301xXsPAg4s+doTuL6X5+SBVOeD4nExBxP0uJ4ETgeeLxZ1dOr8vZAr+08gEgIlKJHq8SJQCwxx98XuvhhYClxOkKR6soygGAIAMzse+Ly7z3X3GQQ9rw3AMSXPGV3yXJHYaIhPpEq4u5vZvcBvzOwsYCUwm6C44qVenvZX4ICSxxngcjNbTlAxeBgwtvhzpwPYdMifSGzUgxKpLp8lSCa/IzjldDjwUXdf18v99xNU+wHg7rcD5xP0ul4Gvgt82d0fKXnOFOABd9cQn8RKC3VFBjAzayA46vwT7v5sGfengX8SVAo+FnJ4In1SD0pkAHP3ZoLe0lllPuU44FUlJ0kCJSiRge9HwP7WrTa9u2Lv6VvAFyKJSmQLNMQnIiKJpB6UiIgkkhKUiIgkkhKUiIgkkhKUiIgkkhKUiIgk0v8Hrt6KFYEAlJ8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "plot_velocity(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Suppose we drop a quarter from the Empire State Building and find that its flight time is 19.1 seconds.  Use this measurement to estimate the terminal velocity.\n",
    "\n",
    "1. You can get the relevant dimensions of a quarter from https://en.wikipedia.org/wiki/Quarter_(United_States_coin).\n",
    "\n",
    "2. Create a `Params` object with the system parameters.  We don't know `v_term`, so we'll start with the inital guess `v_term = 18 * m / s`.\n",
    "\n",
    "3. Use `make_system` to create a `System` object.\n",
    "\n",
    "4. Call `run_ode_solver` to simulate the system.  How does the flight time of the simulation compare to the measurement?\n",
    "\n",
    "5. Try a few different values of `t_term` and see if you can get the simulated flight time close to 19.1 seconds.\n",
    "\n",
    "6. Optionally, write an error function and use `root_scalar` to improve your estimate.\n",
    "\n",
    "7. Use your best estimate of `v_term` to compute `C_d`.\n",
    "\n",
    "Note: I fabricated the observed flight time, so don't take the results of this exercise too seriously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>height</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.00567 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.02426 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>18.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>flight_time</th>\n",
       "      <td>19.1 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "height                         381 meter\n",
       "v_init                0.0 meter / second\n",
       "g                9.8 meter / second ** 2\n",
       "mass                    0.00567 kilogram\n",
       "diameter                   0.02426 meter\n",
       "rho            1.2 kilogram / meter ** 3\n",
       "v_term               18.0 meter / second\n",
       "flight_time                  19.1 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Here's a `Params` object with the dimensions of a quarter,\n",
    "# the observed flight time and our initial guess for `v_term`\n",
    "\n",
    "params3 = Params(params,\n",
    "                 mass = 5.67e-3 * kg,\n",
    "                 diameter = 24.26e-3 * m,\n",
    "                 v_term = 18 * m / s,\n",
    "                 flight_time = 19.1 * s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>height</th>\n",
       "      <td>381 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_init</th>\n",
       "      <td>0.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>g</th>\n",
       "      <td>9.8 meter / second ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mass</th>\n",
       "      <td>0.00567 kilogram</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>diameter</th>\n",
       "      <td>0.02426 meter</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rho</th>\n",
       "      <td>1.2 kilogram / meter ** 3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>v_term</th>\n",
       "      <td>18.0 meter / second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>flight_time</th>\n",
       "      <td>19.1 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>area</th>\n",
       "      <td>0.000462244204111976 meter ** 2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C_d</th>\n",
       "      <td>0.6183600157463346 dimensionless</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>init</th>\n",
       "      <td>y             381 meter\n",
       "v    0.0 meter / secon...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>t_end</th>\n",
       "      <td>30 second</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>dt</th>\n",
       "      <td>0.3 second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "height                                                 381 meter\n",
       "v_init                                        0.0 meter / second\n",
       "g                                        9.8 meter / second ** 2\n",
       "mass                                            0.00567 kilogram\n",
       "diameter                                           0.02426 meter\n",
       "rho                                    1.2 kilogram / meter ** 3\n",
       "v_term                                       18.0 meter / second\n",
       "flight_time                                          19.1 second\n",
       "area                             0.000462244204111976 meter ** 2\n",
       "C_d                             0.6183600157463346 dimensionless\n",
       "init           y             381 meter\n",
       "v    0.0 meter / secon...\n",
       "t_end                                                  30 second\n",
       "dt                                                    0.3 second\n",
       "dtype: object"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Now we can make a `System` object\n",
    "\n",
    "system3 = make_system(params3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>success</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>message</th>\n",
       "      <td>A termination event occurred.</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "success                             True\n",
       "message    A termination event occurred.\n",
       "dtype: object"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Run the simulation\n",
    "\n",
    "results, details = run_ode_solver(system3, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "22.438163885803736 second"
      ],
      "text/latex": [
       "$22.438163885803736\\ \\mathrm{second}$"
      ],
      "text/plain": [
       "22.438163885803736 <Unit('second')>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# And get the flight time\n",
    "\n",
    "flight_time = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution\n",
    "\n",
    "# The flight time is a little long, so we could increase `v_term` and try again.\n",
    "\n",
    "# Or we could write an error function\n",
    "\n",
    "def error_func(guess, params):\n",
    "    \"\"\"Final height as a function of C_d.\n",
    "    \n",
    "    guess: guess at v_term\n",
    "    params: Params object\n",
    "    \n",
    "    returns: height in m\n",
    "    \"\"\"\n",
    "    print(guess)\n",
    "    params = Params(params, v_term=guess)\n",
    "    system = make_system(params)\n",
    "    results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "    flight_time = get_last_label(results) * s\n",
    "    error = flight_time - params.flight_time\n",
    "    return magnitude(error)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18.0 meter / second\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "3.338163885803734"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# We can test the error function like this\n",
    "v_guess1 = 18 * m / s\n",
    "error_func(v_guess1, params3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22.0 meter / second\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "-0.22714064382192944"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "v_guess2 = 22 * m / s\n",
    "error_func(v_guess2, params3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18.0 meter / second\n",
      "22.0 meter / second\n",
      "20.0 meter / second\n",
      "21.0 meter / second\n",
      "21.5 meter / second\n",
      "21.75 meter / second\n",
      "21.625 meter / second\n",
      "21.6875 meter / second\n",
      "21.71875 meter / second\n",
      "21.703125 meter / second\n",
      "21.6953125 meter / second\n",
      "21.69140625 meter / second\n",
      "21.689453125 meter / second\n",
      "21.6884765625 meter / second\n",
      "21.68798828125 meter / second\n",
      "21.687744140625 meter / second\n",
      "21.6878662109375 meter / second\n",
      "21.68792724609375 meter / second\n",
      "21.687896728515625 meter / second\n",
      "21.687881469726562 meter / second\n",
      "21.68787384033203 meter / second\n",
      "21.687877655029297 meter / second\n",
      "21.68787956237793 meter / second\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>values</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>converged</th>\n",
       "      <td>True</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>root</th>\n",
       "      <td>21.687878608703613 meter / second</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "converged                                 True\n",
       "root         21.687878608703613 meter / second\n",
       "dtype: object"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Now we can use `root_scalar` to find the value of `v_term` that yields the measured flight time.\n",
    "\n",
    "res = root_bisect(error_func, [v_guess1, v_guess2], params3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "21.687878608703613 meter/second"
      ],
      "text/latex": [
       "$21.687878608703613\\ \\frac{\\mathrm{meter}}{\\mathrm{second}}$"
      ],
      "text/plain": [
       "21.687878608703613 <Unit('meter / second')>"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "v_term_solution = res.root"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "0.4259437619496639 dimensionless"
      ],
      "text/latex": [
       "$0.4259437619496639\\ dimensionless$"
      ],
      "text/plain": [
       "0.4259437619496639 <Unit('dimensionless')>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Solution\n",
    "\n",
    "# Plugging in the estimated value, we can use `make_system` to compute `C_d`\n",
    "\n",
    "params_solution = Params(params3, v_term=v_term_solution)\n",
    "system = make_system(params_solution)\n",
    "system.C_d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
